GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EIG.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2017 John W. Eaton
4 Copyright (C) 2016 Barbara Lócsi
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #if defined (HAVE_CONFIG_H)
25 # include "config.h"
26 #endif
27 
28 #include "EIG.h"
29 #include "dColVector.h"
30 #include "lo-error.h"
31 #include "lo-lapack-proto.h"
32 
33 octave_idx_type
34 EIG::init (const Matrix& a, bool calc_rev, bool calc_lev, bool balance)
35 {
36  if (a.any_element_is_inf_or_nan ())
37  (*current_liboctave_error_handler)
38  ("EIG: matrix contains Inf or NaN values");
39 
40  if (a.is_symmetric ())
41  return symmetric_init (a, calc_rev, calc_lev);
42 
43  octave_idx_type n = a.rows ();
44 
45  if (n != a.cols ())
46  (*current_liboctave_error_handler) ("EIG requires square matrix");
47 
48  octave_idx_type info = 0;
49 
50  Matrix atmp = a;
51  double *tmp_data = atmp.fortran_vec ();
52 
53  Array<double> wr (dim_vector (n, 1));
54  double *pwr = wr.fortran_vec ();
55 
56  Array<double> wi (dim_vector (n, 1));
57  double *pwi = wi.fortran_vec ();
58 
59  octave_idx_type tnvr = calc_rev ? n : 0;
60  Matrix vr (tnvr, tnvr);
61  double *pvr = vr.fortran_vec ();
62 
63  octave_idx_type tnvl = calc_lev ? n : 0;
64  Matrix vl (tnvl, tnvl);
65  double *pvl = vl.fortran_vec ();
66 
67  octave_idx_type lwork = -1;
68  double dummy_work;
69 
70  octave_idx_type ilo;
71  octave_idx_type ihi;
72 
73  Array<double> scale (dim_vector (n, 1));
74  double *pscale = scale.fortran_vec ();
75 
76  double abnrm;
77 
78  Array<double> rconde (dim_vector (n, 1));
79  double *prconde = rconde.fortran_vec ();
80 
81  Array<double> rcondv (dim_vector (n, 1));
82  double *prcondv = rcondv.fortran_vec ();
83 
84  octave_idx_type dummy_iwork;
85 
86  F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
87  F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
88  F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
89  F77_CONST_CHAR_ARG2 ("N", 1),
90  n, tmp_data, n, pwr, pwi, pvl,
91  n, pvr, n, ilo, ihi, pscale,
92  abnrm, prconde, prcondv, &dummy_work,
93  lwork, &dummy_iwork, info
94  F77_CHAR_ARG_LEN (1)
95  F77_CHAR_ARG_LEN (1)
96  F77_CHAR_ARG_LEN (1)
97  F77_CHAR_ARG_LEN (1)));
98 
99  if (info != 0)
100  (*current_liboctave_error_handler) ("dgeevx workspace query failed");
101 
102  lwork = static_cast<octave_idx_type> (dummy_work);
103  Array<double> work (dim_vector (lwork, 1));
104  double *pwork = work.fortran_vec ();
105 
106  F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
107  F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
108  F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
109  F77_CONST_CHAR_ARG2 ("N", 1),
110  n, tmp_data, n, pwr, pwi, pvl,
111  n, pvr, n, ilo, ihi, pscale,
112  abnrm, prconde, prcondv, pwork,
113  lwork, &dummy_iwork, info
114  F77_CHAR_ARG_LEN (1)
115  F77_CHAR_ARG_LEN (1)
116  F77_CHAR_ARG_LEN (1)
117  F77_CHAR_ARG_LEN (1)));
118 
119  if (info < 0)
120  (*current_liboctave_error_handler) ("unrecoverable error in dgeevx");
121 
122  if (info > 0)
123  (*current_liboctave_error_handler) ("dgeevx failed to converge");
124 
125  lambda.resize (n);
126  octave_idx_type nvr = calc_rev ? n : 0;
127  v.resize (nvr, nvr);
128  octave_idx_type nvl = calc_lev ? n : 0;
129  w.resize (nvl, nvl);
130 
131  for (octave_idx_type j = 0; j < n; j++)
132  {
133  if (wi.elem (j) == 0.0)
134  {
135  lambda.elem (j) = Complex (wr.elem (j));
136  for (octave_idx_type i = 0; i < nvr; i++)
137  v.elem (i, j) = vr.elem (i, j);
138 
139  for (octave_idx_type i = 0; i < nvl; i++)
140  w.elem (i, j) = vl.elem (i, j);
141  }
142  else
143  {
144  if (j+1 >= n)
145  (*current_liboctave_error_handler) ("EIG: internal error");
146 
147  lambda.elem (j) = Complex (wr.elem (j), wi.elem (j));
148  lambda.elem (j+1) = Complex (wr.elem (j+1), wi.elem (j+1));
149 
150  for (octave_idx_type i = 0; i < nvr; i++)
151  {
152  double real_part = vr.elem (i, j);
153  double imag_part = vr.elem (i, j+1);
154  v.elem (i, j) = Complex (real_part, imag_part);
155  v.elem (i, j+1) = Complex (real_part, -imag_part);
156  }
157 
158  for (octave_idx_type i = 0; i < nvl; i++)
159  {
160  double real_part = vl.elem (i, j);
161  double imag_part = vl.elem (i, j+1);
162  w.elem (i, j) = Complex (real_part, imag_part);
163  w.elem (i, j+1) = Complex (real_part, -imag_part);
164  }
165  j++;
166  }
167  }
168 
169  return info;
170 }
171 
172 octave_idx_type
173 EIG::symmetric_init (const Matrix& a, bool calc_rev, bool calc_lev)
174 {
175  octave_idx_type n = a.rows ();
176 
177  if (n != a.cols ())
178  (*current_liboctave_error_handler) ("EIG requires square matrix");
179 
180  octave_idx_type info = 0;
181 
182  Matrix atmp = a;
183  double *tmp_data = atmp.fortran_vec ();
184 
185  ColumnVector wr (n);
186  double *pwr = wr.fortran_vec ();
187 
188  octave_idx_type lwork = -1;
189  double dummy_work;
190 
191  F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
192  F77_CONST_CHAR_ARG2 ("U", 1),
193  n, tmp_data, n, pwr, &dummy_work, lwork, info
194  F77_CHAR_ARG_LEN (1)
195  F77_CHAR_ARG_LEN (1)));
196 
197  if (info != 0)
198  (*current_liboctave_error_handler) ("dsyev workspace query failed");
199 
200  lwork = static_cast<octave_idx_type> (dummy_work);
201  Array<double> work (dim_vector (lwork, 1));
202  double *pwork = work.fortran_vec ();
203 
204  F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
205  F77_CONST_CHAR_ARG2 ("U", 1),
206  n, tmp_data, n, pwr, pwork, lwork, info
207  F77_CHAR_ARG_LEN (1)
208  F77_CHAR_ARG_LEN (1)));
209 
210  if (info < 0)
211  (*current_liboctave_error_handler) ("unrecoverable error in dsyev");
212 
213  if (info > 0)
214  (*current_liboctave_error_handler) ("dsyev failed to converge");
215 
216  lambda = ComplexColumnVector (wr);
217  v = calc_rev ? ComplexMatrix (atmp) : ComplexMatrix ();
218  w = calc_lev ? ComplexMatrix (atmp) : ComplexMatrix ();
219 
220  return info;
221 }
222 
223 octave_idx_type
224 EIG::init (const ComplexMatrix& a, bool calc_rev, bool calc_lev, bool balance)
225 {
226  if (a.any_element_is_inf_or_nan ())
227  (*current_liboctave_error_handler)
228  ("EIG: matrix contains Inf or NaN values");
229 
230  if (a.is_hermitian ())
231  return hermitian_init (a, calc_rev, calc_lev);
232 
233  octave_idx_type n = a.rows ();
234 
235  if (n != a.cols ())
236  (*current_liboctave_error_handler) ("EIG requires square matrix");
237 
238  octave_idx_type info = 0;
239 
240  ComplexMatrix atmp = a;
241  Complex *tmp_data = atmp.fortran_vec ();
242 
243  ComplexColumnVector wr (n);
244  Complex *pw = wr.fortran_vec ();
245 
246  octave_idx_type nvr = calc_rev ? n : 0;
247  ComplexMatrix vrtmp (nvr, nvr);
248  Complex *pvr = vrtmp.fortran_vec ();
249 
250  octave_idx_type nvl = calc_lev ? n : 0;
251  ComplexMatrix vltmp (nvl, nvl);
252  Complex *pvl = vltmp.fortran_vec ();
253 
254  octave_idx_type lwork = -1;
255  Complex dummy_work;
256 
257  octave_idx_type lrwork = 2*n;
258  Array<double> rwork (dim_vector (lrwork, 1));
259  double *prwork = rwork.fortran_vec ();
260 
261  octave_idx_type ilo;
262  octave_idx_type ihi;
263 
264  Array<double> scale (dim_vector (n, 1));
265  double *pscale = scale.fortran_vec ();
266 
267  double abnrm;
268 
269  Array<double> rconde (dim_vector (n, 1));
270  double *prconde = rconde.fortran_vec ();
271 
272  Array<double> rcondv (dim_vector (n, 1));
273  double *prcondv = rcondv.fortran_vec ();
274 
275  F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
276  F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
277  F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
278  F77_CONST_CHAR_ARG2 ("N", 1),
279  n, F77_DBLE_CMPLX_ARG (tmp_data), n,
280  F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (pvl), n,
281  F77_DBLE_CMPLX_ARG (pvr), n, ilo, ihi, pscale, abnrm,
282  prconde, prcondv,
283  F77_DBLE_CMPLX_ARG (&dummy_work), lwork, prwork, info
284  F77_CHAR_ARG_LEN (1)
285  F77_CHAR_ARG_LEN (1)
286  F77_CHAR_ARG_LEN (1)
287  F77_CHAR_ARG_LEN (1)));
288 
289  if (info != 0)
290  (*current_liboctave_error_handler) ("zgeevx workspace query failed");
291 
292  lwork = static_cast<octave_idx_type> (dummy_work.real ());
293  Array<Complex> work (dim_vector (lwork, 1));
294  Complex *pwork = work.fortran_vec ();
295 
296  F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
297  F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
298  F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
299  F77_CONST_CHAR_ARG2 ("N", 1),
300  n, F77_DBLE_CMPLX_ARG (tmp_data), n,
301  F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (pvl), n,
302  F77_DBLE_CMPLX_ARG (pvr), n, ilo, ihi, pscale, abnrm,
303  prconde, prcondv,
304  F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info
305  F77_CHAR_ARG_LEN (1)
306  F77_CHAR_ARG_LEN (1)
307  F77_CHAR_ARG_LEN (1)
308  F77_CHAR_ARG_LEN (1)));
309 
310  if (info < 0)
311  (*current_liboctave_error_handler) ("unrecoverable error in zgeevx");
312 
313  if (info > 0)
314  (*current_liboctave_error_handler) ("zgeevx failed to converge");
315 
316  lambda = wr;
317  v = vrtmp;
318  w = vltmp;
319 
320  return info;
321 }
322 
323 octave_idx_type
324 EIG::hermitian_init (const ComplexMatrix& a, bool calc_rev, bool calc_lev)
325 {
326  octave_idx_type n = a.rows ();
327 
328  if (n != a.cols ())
329  (*current_liboctave_error_handler) ("EIG requires square matrix");
330 
331  octave_idx_type info = 0;
332 
333  ComplexMatrix atmp = a;
334  Complex *tmp_data = atmp.fortran_vec ();
335 
336  ColumnVector wr (n);
337  double *pwr = wr.fortran_vec ();
338 
339  octave_idx_type lwork = -1;
340  Complex dummy_work;
341 
342  octave_idx_type lrwork = 3*n;
343  Array<double> rwork (dim_vector (lrwork, 1));
344  double *prwork = rwork.fortran_vec ();
345 
346  F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
347  F77_CONST_CHAR_ARG2 ("U", 1),
348  n, F77_DBLE_CMPLX_ARG (tmp_data), n, pwr,
349  F77_DBLE_CMPLX_ARG (&dummy_work), lwork,
350  prwork, info
351  F77_CHAR_ARG_LEN (1)
352  F77_CHAR_ARG_LEN (1)));
353 
354  if (info != 0)
355  (*current_liboctave_error_handler) ("zheev workspace query failed");
356 
357  lwork = static_cast<octave_idx_type> (dummy_work.real ());
358  Array<Complex> work (dim_vector (lwork, 1));
359  Complex *pwork = work.fortran_vec ();
360 
361  F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
362  F77_CONST_CHAR_ARG2 ("U", 1),
363  n, F77_DBLE_CMPLX_ARG (tmp_data), n, pwr,
364  F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info
365  F77_CHAR_ARG_LEN (1)
366  F77_CHAR_ARG_LEN (1)));
367 
368  if (info < 0)
369  (*current_liboctave_error_handler) ("unrecoverable error in zheev");
370 
371  if (info > 0)
372  (*current_liboctave_error_handler) ("zheev failed to converge");
373 
374  lambda = ComplexColumnVector (wr);
375  v = calc_rev ? ComplexMatrix (atmp) : ComplexMatrix ();
376  w = calc_lev ? ComplexMatrix (atmp) : ComplexMatrix ();
377 
378  return info;
379 }
380 
381 octave_idx_type
382 EIG::init (const Matrix& a, const Matrix& b, bool calc_rev, bool calc_lev,
383  bool force_qz)
384 {
385  if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
386  (*current_liboctave_error_handler)
387  ("EIG: matrix contains Inf or NaN values");
388 
389  octave_idx_type n = a.rows ();
390  octave_idx_type nb = b.rows ();
391 
392  if (n != a.cols () || nb != b.cols ())
393  (*current_liboctave_error_handler) ("EIG requires square matrix");
394 
395  if (n != nb)
396  (*current_liboctave_error_handler) ("EIG requires same size matrices");
397 
398  octave_idx_type info = 0;
399 
400  Matrix tmp = b;
401  double *tmp_data = tmp.fortran_vec ();
402 
403  if (! force_qz)
404  {
405  F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
406  n, tmp_data, n,
407  info
408  F77_CHAR_ARG_LEN (1)));
409 
410  if (a.is_symmetric () && b.is_symmetric () && info == 0)
411  return symmetric_init (a, b, calc_rev, calc_lev);
412  }
413 
414  Matrix atmp = a;
415  double *atmp_data = atmp.fortran_vec ();
416 
417  Matrix btmp = b;
418  double *btmp_data = btmp.fortran_vec ();
419 
420  Array<double> ar (dim_vector (n, 1));
421  double *par = ar.fortran_vec ();
422 
423  Array<double> ai (dim_vector (n, 1));
424  double *pai = ai.fortran_vec ();
425 
426  Array<double> beta (dim_vector (n, 1));
427  double *pbeta = beta.fortran_vec ();
428 
429  octave_idx_type tnvr = calc_rev ? n : 0;
430  Matrix vr (tnvr, tnvr);
431  double *pvr = vr.fortran_vec ();
432 
433  octave_idx_type tnvl = calc_lev ? n : 0;
434  Matrix vl (tnvl, tnvl);
435  double *pvl = vl.fortran_vec ();
436 
437  octave_idx_type lwork = -1;
438  double dummy_work;
439 
440  F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
441  F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
442  n, atmp_data, n, btmp_data, n,
443  par, pai, pbeta,
444  pvl, n, pvr, n,
445  &dummy_work, lwork, info
446  F77_CHAR_ARG_LEN (1)
447  F77_CHAR_ARG_LEN (1)));
448 
449  if (info != 0)
450  (*current_liboctave_error_handler) ("dggev workspace query failed");
451 
452  lwork = static_cast<octave_idx_type> (dummy_work);
453  Array<double> work (dim_vector (lwork, 1));
454  double *pwork = work.fortran_vec ();
455 
456  F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
457  F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
458  n, atmp_data, n, btmp_data, n,
459  par, pai, pbeta,
460  pvl, n, pvr, n,
461  pwork, lwork, info
462  F77_CHAR_ARG_LEN (1)
463  F77_CHAR_ARG_LEN (1)));
464 
465  if (info < 0)
466  (*current_liboctave_error_handler) ("unrecoverable error in dggev");
467 
468  if (info > 0)
469  (*current_liboctave_error_handler) ("dggev failed to converge");
470 
471  lambda.resize (n);
472  octave_idx_type nvr = calc_rev ? n : 0;
473  v.resize (nvr, nvr);
474 
475  octave_idx_type nvl = calc_lev ? n : 0;
476  w.resize (nvl, nvl);
477 
478  for (octave_idx_type j = 0; j < n; j++)
479  {
480  if (ai.elem (j) == 0.0)
481  {
482  lambda.elem (j) = Complex (ar.elem (j) / beta.elem (j));
483  for (octave_idx_type i = 0; i < nvr; i++)
484  v.elem (i, j) = vr.elem (i, j);
485  for (octave_idx_type i = 0; i < nvl; i++)
486  w.elem (i, j) = vl.elem (i, j);
487  }
488  else
489  {
490  if (j+1 >= n)
491  (*current_liboctave_error_handler) ("EIG: internal error");
492 
493  lambda.elem (j) = Complex (ar.elem (j) / beta.elem (j),
494  ai.elem (j) / beta.elem (j));
495  lambda.elem (j+1) = Complex (ar.elem (j+1) / beta.elem (j+1),
496  ai.elem (j+1) / beta.elem (j+1));
497 
498  for (octave_idx_type i = 0; i < nvr; i++)
499  {
500  double real_part = vr.elem (i, j);
501  double imag_part = vr.elem (i, j+1);
502  v.elem (i, j) = Complex (real_part, imag_part);
503  v.elem (i, j+1) = Complex (real_part, -imag_part);
504  }
505  for (octave_idx_type i = 0; i < nvl; i++)
506  {
507  double real_part = vl.elem (i, j);
508  double imag_part = vl.elem (i, j+1);
509  w.elem (i, j) = Complex (real_part, imag_part);
510  w.elem (i, j+1) = Complex (real_part, -imag_part);
511  }
512  j++;
513  }
514  }
515 
516  return info;
517 }
518 
519 octave_idx_type
520 EIG::symmetric_init (const Matrix& a, const Matrix& b, bool calc_rev,
521  bool calc_lev)
522 {
523  octave_idx_type n = a.rows ();
524  octave_idx_type nb = b.rows ();
525 
526  if (n != a.cols () || nb != b.cols ())
527  (*current_liboctave_error_handler) ("EIG requires square matrix");
528 
529  if (n != nb)
530  (*current_liboctave_error_handler) ("EIG requires same size matrices");
531 
532  octave_idx_type info = 0;
533 
534  Matrix atmp = a;
535  double *atmp_data = atmp.fortran_vec ();
536 
537  Matrix btmp = b;
538  double *btmp_data = btmp.fortran_vec ();
539 
540  ColumnVector wr (n);
541  double *pwr = wr.fortran_vec ();
542 
543  octave_idx_type lwork = -1;
544  double dummy_work;
545 
546  F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
547  F77_CONST_CHAR_ARG2 ("U", 1),
548  n, atmp_data, n,
549  btmp_data, n,
550  pwr, &dummy_work, lwork, info
551  F77_CHAR_ARG_LEN (1)
552  F77_CHAR_ARG_LEN (1)));
553 
554  if (info != 0)
555  (*current_liboctave_error_handler) ("dsygv workspace query failed");
556 
557  lwork = static_cast<octave_idx_type> (dummy_work);
558  Array<double> work (dim_vector (lwork, 1));
559  double *pwork = work.fortran_vec ();
560 
561  F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
562  F77_CONST_CHAR_ARG2 ("U", 1),
563  n, atmp_data, n,
564  btmp_data, n,
565  pwr, pwork, lwork, info
566  F77_CHAR_ARG_LEN (1)
567  F77_CHAR_ARG_LEN (1)));
568 
569  if (info < 0)
570  (*current_liboctave_error_handler) ("unrecoverable error in dsygv");
571 
572  if (info > 0)
573  (*current_liboctave_error_handler) ("dsygv failed to converge");
574 
575  lambda = ComplexColumnVector (wr);
576  v = calc_rev ? ComplexMatrix (atmp) : ComplexMatrix ();
577  w = calc_lev ? ComplexMatrix (atmp) : ComplexMatrix ();
578 
579  return info;
580 }
581 
582 octave_idx_type
583 EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_rev,
584  bool calc_lev, bool force_qz)
585 {
586  if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
587  (*current_liboctave_error_handler)
588  ("EIG: matrix contains Inf or NaN values");
589 
590  octave_idx_type n = a.rows ();
591  octave_idx_type nb = b.rows ();
592 
593  if (n != a.cols () || nb != b.cols ())
594  (*current_liboctave_error_handler) ("EIG requires square matrix");
595 
596  if (n != nb)
597  (*current_liboctave_error_handler) ("EIG requires same size matrices");
598 
599  octave_idx_type info = 0;
600 
601  ComplexMatrix tmp = b;
602  Complex*tmp_data = tmp.fortran_vec ();
603 
604  if (! force_qz)
605  {
606  F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
607  n, F77_DBLE_CMPLX_ARG (tmp_data), n,
608  info
609  F77_CHAR_ARG_LEN (1)));
610 
611  if (a.is_hermitian () && b.is_hermitian () && info == 0)
612  return hermitian_init (a, b, calc_rev, calc_lev);
613  }
614 
615  ComplexMatrix atmp = a;
616  Complex *atmp_data = atmp.fortran_vec ();
617 
618  ComplexMatrix btmp = b;
619  Complex *btmp_data = btmp.fortran_vec ();
620 
621  ComplexColumnVector alpha (n);
622  Complex *palpha = alpha.fortran_vec ();
623 
624  ComplexColumnVector beta (n);
625  Complex *pbeta = beta.fortran_vec ();
626 
627  octave_idx_type nvr = calc_rev ? n : 0;
628  ComplexMatrix vrtmp (nvr, nvr);
629  Complex *pvr = vrtmp.fortran_vec ();
630 
631  octave_idx_type nvl = calc_lev ? n : 0;
632  ComplexMatrix vltmp (nvl, nvl);
633  Complex *pvl = vltmp.fortran_vec ();
634 
635  octave_idx_type lwork = -1;
636  Complex dummy_work;
637 
638  octave_idx_type lrwork = 8*n;
639  Array<double> rwork (dim_vector (lrwork, 1));
640  double *prwork = rwork.fortran_vec ();
641 
642  F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
643  F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
644  n, F77_DBLE_CMPLX_ARG (atmp_data), n,
645  F77_DBLE_CMPLX_ARG (btmp_data), n,
646  F77_DBLE_CMPLX_ARG (palpha),
647  F77_DBLE_CMPLX_ARG (pbeta),
648  F77_DBLE_CMPLX_ARG (pvl), n,
649  F77_DBLE_CMPLX_ARG (pvr), n,
650  F77_DBLE_CMPLX_ARG (&dummy_work), lwork, prwork, info
651  F77_CHAR_ARG_LEN (1)
652  F77_CHAR_ARG_LEN (1)));
653 
654  if (info != 0)
655  (*current_liboctave_error_handler) ("zggev workspace query failed");
656 
657  lwork = static_cast<octave_idx_type> (dummy_work.real ());
658  Array<Complex> work (dim_vector (lwork, 1));
659  Complex *pwork = work.fortran_vec ();
660 
661  F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
662  F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
663  n, F77_DBLE_CMPLX_ARG (atmp_data), n,
664  F77_DBLE_CMPLX_ARG (btmp_data), n,
665  F77_DBLE_CMPLX_ARG (palpha),
666  F77_DBLE_CMPLX_ARG (pbeta),
667  F77_DBLE_CMPLX_ARG (pvl), n,
668  F77_DBLE_CMPLX_ARG (pvr), n,
669  F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info
670  F77_CHAR_ARG_LEN (1)
671  F77_CHAR_ARG_LEN (1)));
672 
673  if (info < 0)
674  (*current_liboctave_error_handler) ("unrecoverable error in zggev");
675 
676  if (info > 0)
677  (*current_liboctave_error_handler) ("zggev failed to converge");
678 
679  lambda.resize (n);
680 
681  for (octave_idx_type j = 0; j < n; j++)
682  lambda.elem (j) = alpha.elem (j) / beta.elem (j);
683 
684  v = vrtmp;
685  w = vltmp;
686 
687  return info;
688 }
689 
690 octave_idx_type
691 EIG::hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b,
692  bool calc_rev, bool calc_lev)
693 {
694  octave_idx_type n = a.rows ();
695  octave_idx_type nb = b.rows ();
696 
697  if (n != a.cols () || nb != b.cols ())
698  (*current_liboctave_error_handler) ("EIG requires square matrix");
699 
700  if (n != nb)
701  (*current_liboctave_error_handler) ("EIG requires same size matrices");
702 
703  octave_idx_type info = 0;
704 
705  ComplexMatrix atmp = a;
706  Complex *atmp_data = atmp.fortran_vec ();
707 
708  ComplexMatrix btmp = b;
709  Complex *btmp_data = btmp.fortran_vec ();
710 
711  ColumnVector wr (n);
712  double *pwr = wr.fortran_vec ();
713 
714  octave_idx_type lwork = -1;
715  Complex dummy_work;
716 
717  octave_idx_type lrwork = 3*n;
718  Array<double> rwork (dim_vector (lrwork, 1));
719  double *prwork = rwork.fortran_vec ();
720 
721  F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
722  F77_CONST_CHAR_ARG2 ("U", 1),
723  n, F77_DBLE_CMPLX_ARG (atmp_data), n,
724  F77_DBLE_CMPLX_ARG (btmp_data), n,
725  pwr, F77_DBLE_CMPLX_ARG (&dummy_work), lwork,
726  prwork, info
727  F77_CHAR_ARG_LEN (1)
728  F77_CHAR_ARG_LEN (1)));
729 
730  if (info != 0)
731  (*current_liboctave_error_handler) ("zhegv workspace query failed");
732 
733  lwork = static_cast<octave_idx_type> (dummy_work.real ());
734  Array<Complex> work (dim_vector (lwork, 1));
735  Complex *pwork = work.fortran_vec ();
736 
737  F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
738  F77_CONST_CHAR_ARG2 ("U", 1),
739  n, F77_DBLE_CMPLX_ARG (atmp_data), n,
740  F77_DBLE_CMPLX_ARG (btmp_data), n,
741  pwr, F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info
742  F77_CHAR_ARG_LEN (1)
743  F77_CHAR_ARG_LEN (1)));
744 
745  if (info < 0)
746  (*current_liboctave_error_handler) ("unrecoverable error in zhegv");
747 
748  if (info > 0)
749  (*current_liboctave_error_handler) ("zhegv failed to converge");
750 
751  lambda = ComplexColumnVector (wr);
752  v = calc_rev ? ComplexMatrix (atmp) : ComplexMatrix ();
753  w = calc_lev ? ComplexMatrix (atmp) : ComplexMatrix ();
754 
755  return info;
756 }
Matrix::is_symmetric
bool is_symmetric(void) const
Definition: dMatrix.cc:137
ComplexColumnVector::resize
void resize(octave_idx_type n, const Complex &rfv=Complex(0))
Definition: CColVector.h:141
current_liboctave_error_handler
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:345
Array::elem
T & elem(octave_idx_type n)
Definition: Array.h:482
F77_XFCN
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:52
EIG::symmetric_init
octave_idx_type symmetric_init(const Matrix &a, bool calc_rev, bool calc_lev)
Definition: EIG.cc:173
Array::rows
octave_idx_type rows(void) const
Definition: Array.h:401
a
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
EIG::hermitian_init
octave_idx_type hermitian_init(const ComplexMatrix &a, bool calc_rev, bool calc_lev)
Definition: EIG.cc:324
ComplexMatrix::is_hermitian
bool is_hermitian(void) const
Definition: CMatrix.cc:177
EIG::init
octave_idx_type init(const Matrix &a, bool calc_rev, bool calc_lev, bool balance)
Definition: EIG.cc:34
tmp
double tmp
Definition: data.cc:6300
Matrix
Definition: dMatrix.h:37
ComplexNDArray::any_element_is_inf_or_nan
bool any_element_is_inf_or_nan(void) const
Definition: CNDArray.cc:508
EIG::w
ComplexMatrix w
Definition: EIG.h:127
EIG::lambda
ComplexColumnVector lambda
Definition: EIG.h:125
EIG::ComplexMatrix
friend class ComplexMatrix
Definition: EIG.h:39
ComplexMatrix::resize
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
Definition: CMatrix.h:184
wi
static double wi[256]
Definition: randmtzig.cc:440
i
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
EIG::v
ComplexMatrix v
Definition: EIG.h:126
b
b
Definition: cellfun.cc:398
scale
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5150
NDArray::any_element_is_inf_or_nan
bool any_element_is_inf_or_nan(void) const
Definition: dNDArray.cc:561
octave_idx_type
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:31
Array::fortran_vec
const T * fortran_vec(void) const
Definition: Array.h:584
ar
ar
Definition: __qp__.cc:492
Array::cols
octave_idx_type cols(void) const
Definition: Array.h:409
dim_vector
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
Generated on Wed May 10 2017 15:42:57 for GNU Octave by   doxygen 1.8.8