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