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