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