GNU Octave  8.1.0
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
fEIG.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1994-2023 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 
38 FloatEIG::init (const FloatMatrix& a, bool calc_rev, bool calc_lev,
39  bool balance)
40 {
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 
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  m_lambda.resize (n);
132  m_v.resize (nvr, nvr);
133  m_w.resize (nvl, nvl);
134 
135  for (F77_INT j = 0; j < n; j++)
136  {
137  if (wi.elem (j) == 0.0)
138  {
139  m_lambda.elem (j) = FloatComplex (wr.elem (j));
140  for (octave_idx_type i = 0; i < nvr; i++)
141  m_v.elem (i, j) = vr.elem (i, j);
142 
143  for (F77_INT i = 0; i < nvl; i++)
144  m_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  m_lambda.elem (j) = FloatComplex (wr.elem (j), wi.elem (j));
152  m_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  m_v.elem (i, j) = FloatComplex (real_part, imag_part);
159  m_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  m_w.elem (i, j) = FloatComplex (real_part, imag_part);
166  m_w.elem (i, j+1) = FloatComplex (real_part, -imag_part);
167  }
168  j++;
169  }
170  }
171 
172  return info;
173 }
174 
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 
221  m_v = (calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
222  m_w = (calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
223 
224  return info;
225 }
226 
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 
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 
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  m_lambda = wr;
321  m_v = vrtmp;
322  m_w = vltmp;
323 
324  return info;
325 }
326 
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 
381  m_v = (calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
382  m_w = (calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
383 
384  return info;
385 }
386 
388 FloatEIG::init (const FloatMatrix& a, const FloatMatrix& b, bool calc_rev,
389  bool calc_lev, bool force_qz)
390 {
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  m_lambda.resize (n);
480  m_v.resize (nvr, nvr);
481  m_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  m_lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j));
489  for (F77_INT i = 0; i < nvr; i++)
490  m_v.elem (i, j) = vr.elem (i, j);
491 
492  for (F77_INT i = 0; i < nvl; i++)
493  m_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  m_lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j),
501  ai.elem (j) / beta.elem (j));
502  m_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  m_v.elem (i, j) = FloatComplex (real_part, imag_part);
510  m_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  m_w.elem (i, j) = FloatComplex (real_part, imag_part);
517  m_w.elem (i, j+1) = FloatComplex (real_part, -imag_part);
518  }
519  j++;
520  }
521  }
522 
523  return info;
524 }
525 
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 
586  m_v = (calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
587  m_w = (calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
588 
589  return info;
590 }
591 
594  bool calc_rev, bool calc_lev, bool force_qz)
595 {
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 
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  m_lambda.resize (n);
689 
690  for (F77_INT j = 0; j < n; j++)
691  m_lambda.elem (j) = alpha.elem (j) / beta.elem (j);
692 
693  m_v = vrtmp;
694  m_w = vltmp;
695 
696  return info;
697 }
698 
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 
765  m_v = (calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
766  m_w = (calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
767 
768  return info;
769 }
N Dimensional Array with copy-on-write semantics.
Definition: Array.h:130
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type rows(void) const
Definition: Array.h:459
OCTARRAY_API T * fortran_vec(void)
Size of the specified dimension.
Definition: Array-base.cc:1766
OCTARRAY_OVERRIDABLE_FUNC_API T & elem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:562
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type cols(void) const
Definition: Array.h:469
void resize(octave_idx_type n, const FloatComplex &rfv=FloatComplex(0))
Definition: fCColVector.h:152
OCTAVE_API bool ishermitian(void) const
Definition: fCMatrix.cc:174
void resize(octave_idx_type nr, octave_idx_type nc, const FloatComplex &rfv=FloatComplex(0))
Definition: fCMatrix.h:201
OCTAVE_API bool any_element_is_inf_or_nan(void) const
Definition: fCNDArray.cc:271
FloatComplexColumnVector m_lambda
Definition: fEIG.h:129
octave_idx_type symmetric_init(const FloatMatrix &a, bool calc_rev, bool calc_lev)
Definition: fEIG.cc:176
friend class FloatComplexMatrix
Definition: fEIG.h:43
FloatComplexMatrix m_v
Definition: fEIG.h:130
FloatComplexMatrix m_w
Definition: fEIG.h:131
octave_idx_type init(const FloatMatrix &a, bool calc_rev, bool calc_lev, bool balance)
Definition: fEIG.cc:38
octave_idx_type hermitian_init(const FloatComplexMatrix &a, bool calc_rev, bool calc_lev)
Definition: fEIG.cc:328
OCTAVE_API bool issymmetric(void) const
Definition: fMatrix.cc:138
OCTAVE_API bool any_element_is_inf_or_nan(void) const
Definition: fNDArray.cc:281
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:310
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:45
octave_f77_int_type F77_INT
Definition: f77-fcn.h:306
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5851
class OCTAVE_API FloatComplexColumnVector
Definition: mx-fwd.h:48
octave_idx_type n
Definition: mx-inlines.cc:753
std::complex< float > FloatComplex
Definition: oct-cmplx.h:34
static double wi[256]
Definition: randmtzig.cc:464