GNU Octave  3.8.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-2013 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include "fEIG.h"
28 #include "fColVector.h"
29 #include "f77-fcn.h"
30 #include "lo-error.h"
31 
32 extern "C"
33 {
34  F77_RET_T
35  F77_FUNC (sgeev, SGEEV) (F77_CONST_CHAR_ARG_DECL,
37  const octave_idx_type&, float*,
38  const octave_idx_type&, float*, float*, float*,
39  const octave_idx_type&, float*,
40  const octave_idx_type&, float*,
41  const octave_idx_type&, octave_idx_type&
44 
45  F77_RET_T
46  F77_FUNC (cgeev, CGEEV) (F77_CONST_CHAR_ARG_DECL,
48  const octave_idx_type&, FloatComplex*,
49  const octave_idx_type&, FloatComplex*, FloatComplex*,
50  const octave_idx_type&, FloatComplex*,
51  const octave_idx_type&, FloatComplex*,
52  const octave_idx_type&, float*, octave_idx_type&
55 
56  F77_RET_T
57  F77_FUNC (ssyev, SSYEV) (F77_CONST_CHAR_ARG_DECL,
59  const octave_idx_type&, float*,
60  const octave_idx_type&, float*, float*,
61  const octave_idx_type&, octave_idx_type&
64 
65  F77_RET_T
66  F77_FUNC (cheev, CHEEV) (F77_CONST_CHAR_ARG_DECL,
68  const octave_idx_type&, FloatComplex*,
69  const octave_idx_type&, float*, FloatComplex*,
70  const octave_idx_type&, float*, octave_idx_type&
73 
74  F77_RET_T
75  F77_FUNC (spotrf, SPOTRF) (F77_CONST_CHAR_ARG_DECL,
76  const octave_idx_type&, float*,
77  const octave_idx_type&, octave_idx_type&
80 
81  F77_RET_T
82  F77_FUNC (cpotrf, CPOTRF) (F77_CONST_CHAR_ARG_DECL,
83  const octave_idx_type&, FloatComplex*,
84  const octave_idx_type&, octave_idx_type&
87 
88  F77_RET_T
89  F77_FUNC (sggev, SGGEV) (F77_CONST_CHAR_ARG_DECL,
91  const octave_idx_type&, float*,
92  const octave_idx_type&, float*,
93  const octave_idx_type&, float*, float*, float*,
94  float*, const octave_idx_type&, float*,
95  const octave_idx_type&, float*,
96  const octave_idx_type&, octave_idx_type&
99 
100  F77_RET_T
101  F77_FUNC (ssygv, SSYGV) (const octave_idx_type&,
104  const octave_idx_type&, float*,
105  const octave_idx_type&, float*,
106  const octave_idx_type&, float*, float*,
107  const octave_idx_type&, octave_idx_type&
110 
111  F77_RET_T
112  F77_FUNC (cggev, CGGEV) (F77_CONST_CHAR_ARG_DECL,
114  const octave_idx_type&, FloatComplex*,
115  const octave_idx_type&, FloatComplex*,
116  const octave_idx_type&, FloatComplex*,
117  FloatComplex*, FloatComplex*,
118  const octave_idx_type&, FloatComplex*,
119  const octave_idx_type&, FloatComplex*,
120  const octave_idx_type&, float*, octave_idx_type&
123 
124  F77_RET_T
125  F77_FUNC (chegv, CHEGV) (const octave_idx_type&,
128  const octave_idx_type&, FloatComplex*,
129  const octave_idx_type&, FloatComplex*,
130  const octave_idx_type&, float*, FloatComplex*,
131  const octave_idx_type&, float*, octave_idx_type&
134 }
135 
137 FloatEIG::init (const FloatMatrix& a, bool calc_ev)
138 {
139  if (a.any_element_is_inf_or_nan ())
140  {
141  (*current_liboctave_error_handler)
142  ("EIG: matrix contains Inf or NaN values");
143  return -1;
144  }
145 
146  if (a.is_symmetric ())
147  return symmetric_init (a, calc_ev);
148 
149  octave_idx_type n = a.rows ();
150 
151  if (n != a.cols ())
152  {
153  (*current_liboctave_error_handler) ("EIG requires square matrix");
154  return -1;
155  }
156 
157  octave_idx_type info = 0;
158 
159  FloatMatrix atmp = a;
160  float *tmp_data = atmp.fortran_vec ();
161 
162  Array<float> wr (dim_vector (n, 1));
163  float *pwr = wr.fortran_vec ();
164 
165  Array<float> wi (dim_vector (n, 1));
166  float *pwi = wi.fortran_vec ();
167 
168  volatile octave_idx_type nvr = calc_ev ? n : 0;
169  FloatMatrix vr (nvr, nvr);
170  float *pvr = vr.fortran_vec ();
171 
172  octave_idx_type lwork = -1;
173  float dummy_work;
174 
175  float *dummy = 0;
176  octave_idx_type idummy = 1;
177 
178  F77_XFCN (sgeev, SGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
179  F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
180  n, tmp_data, n, pwr, pwi, dummy,
181  idummy, pvr, n, &dummy_work, lwork, info
182  F77_CHAR_ARG_LEN (1)
183  F77_CHAR_ARG_LEN (1)));
184 
185  if (info == 0)
186  {
187  lwork = static_cast<octave_idx_type> (dummy_work);
188  Array<float> work (dim_vector (lwork, 1));
189  float *pwork = work.fortran_vec ();
190 
191  F77_XFCN (sgeev, SGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
192  F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
193  n, tmp_data, n, pwr, pwi, dummy,
194  idummy, pvr, n, pwork, lwork, info
195  F77_CHAR_ARG_LEN (1)
196  F77_CHAR_ARG_LEN (1)));
197 
198  if (info < 0)
199  {
200  (*current_liboctave_error_handler) ("unrecoverable error in sgeev");
201  return info;
202  }
203 
204  if (info > 0)
205  {
206  (*current_liboctave_error_handler) ("sgeev failed to converge");
207  return info;
208  }
209 
210  lambda.resize (n);
211  v.resize (nvr, nvr);
212 
213  for (octave_idx_type j = 0; j < n; j++)
214  {
215  if (wi.elem (j) == 0.0)
216  {
217  lambda.elem (j) = FloatComplex (wr.elem (j));
218  for (octave_idx_type i = 0; i < nvr; i++)
219  v.elem (i, j) = vr.elem (i, j);
220  }
221  else
222  {
223  if (j+1 >= n)
224  {
225  (*current_liboctave_error_handler) ("EIG: internal error");
226  return -1;
227  }
228 
229  lambda.elem (j) = FloatComplex (wr.elem (j), wi.elem (j));
230  lambda.elem (j+1) = FloatComplex (wr.elem (j+1), wi.elem (j+1));
231 
232  for (octave_idx_type i = 0; i < nvr; i++)
233  {
234  float real_part = vr.elem (i, j);
235  float imag_part = vr.elem (i, j+1);
236  v.elem (i, j) = FloatComplex (real_part, imag_part);
237  v.elem (i, j+1) = FloatComplex (real_part, -imag_part);
238  }
239  j++;
240  }
241  }
242  }
243  else
244  (*current_liboctave_error_handler) ("sgeev workspace query failed");
245 
246  return info;
247 }
248 
250 FloatEIG::symmetric_init (const FloatMatrix& a, bool calc_ev)
251 {
252  octave_idx_type n = a.rows ();
253 
254  if (n != a.cols ())
255  {
256  (*current_liboctave_error_handler) ("EIG requires square matrix");
257  return -1;
258  }
259 
260  octave_idx_type info = 0;
261 
262  FloatMatrix atmp = a;
263  float *tmp_data = atmp.fortran_vec ();
264 
265  FloatColumnVector wr (n);
266  float *pwr = wr.fortran_vec ();
267 
268  octave_idx_type lwork = -1;
269  float dummy_work;
270 
271  F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
272  F77_CONST_CHAR_ARG2 ("U", 1),
273  n, tmp_data, n, pwr, &dummy_work, lwork, info
274  F77_CHAR_ARG_LEN (1)
275  F77_CHAR_ARG_LEN (1)));
276 
277  if (info == 0)
278  {
279  lwork = static_cast<octave_idx_type> (dummy_work);
280  Array<float> work (dim_vector (lwork, 1));
281  float *pwork = work.fortran_vec ();
282 
283  F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
284  F77_CONST_CHAR_ARG2 ("U", 1),
285  n, tmp_data, n, pwr, pwork, lwork, info
286  F77_CHAR_ARG_LEN (1)
287  F77_CHAR_ARG_LEN (1)));
288 
289  if (info < 0)
290  {
291  (*current_liboctave_error_handler) ("unrecoverable error in ssyev");
292  return info;
293  }
294 
295  if (info > 0)
296  {
297  (*current_liboctave_error_handler) ("ssyev failed to converge");
298  return info;
299  }
300 
302  v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
303  }
304  else
305  (*current_liboctave_error_handler) ("ssyev workspace query failed");
306 
307  return info;
308 }
309 
311 FloatEIG::init (const FloatComplexMatrix& a, bool calc_ev)
312 {
313  if (a.any_element_is_inf_or_nan ())
314  {
315  (*current_liboctave_error_handler)
316  ("EIG: matrix contains Inf or NaN values");
317  return -1;
318  }
319 
320  if (a.is_hermitian ())
321  return hermitian_init (a, calc_ev);
322 
323  octave_idx_type n = a.rows ();
324 
325  if (n != a.cols ())
326  {
327  (*current_liboctave_error_handler) ("EIG requires square matrix");
328  return -1;
329  }
330 
331  octave_idx_type info = 0;
332 
333  FloatComplexMatrix atmp = a;
334  FloatComplex *tmp_data = atmp.fortran_vec ();
335 
337  FloatComplex *pw = w.fortran_vec ();
338 
339  octave_idx_type nvr = calc_ev ? n : 0;
340  FloatComplexMatrix vtmp (nvr, nvr);
341  FloatComplex *pv = vtmp.fortran_vec ();
342 
343  octave_idx_type lwork = -1;
344  FloatComplex dummy_work;
345 
346  octave_idx_type lrwork = 2*n;
347  Array<float> rwork (dim_vector (lrwork, 1));
348  float *prwork = rwork.fortran_vec ();
349 
350  FloatComplex *dummy = 0;
351  octave_idx_type idummy = 1;
352 
353  F77_XFCN (cgeev, CGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
354  F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
355  n, tmp_data, n, pw, dummy, idummy,
356  pv, n, &dummy_work, lwork, prwork, info
357  F77_CHAR_ARG_LEN (1)
358  F77_CHAR_ARG_LEN (1)));
359 
360  if (info == 0)
361  {
362  lwork = static_cast<octave_idx_type> (dummy_work.real ());
363  Array<FloatComplex> work (dim_vector (lwork, 1));
364  FloatComplex *pwork = work.fortran_vec ();
365 
366  F77_XFCN (cgeev, CGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
367  F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
368  n, tmp_data, n, pw, dummy, idummy,
369  pv, n, pwork, lwork, prwork, info
370  F77_CHAR_ARG_LEN (1)
371  F77_CHAR_ARG_LEN (1)));
372 
373  if (info < 0)
374  {
375  (*current_liboctave_error_handler) ("unrecoverable error in cgeev");
376  return info;
377  }
378 
379  if (info > 0)
380  {
381  (*current_liboctave_error_handler) ("cgeev failed to converge");
382  return info;
383  }
384 
385  lambda = w;
386  v = vtmp;
387  }
388  else
389  (*current_liboctave_error_handler) ("cgeev workspace query failed");
390 
391  return info;
392 }
393 
396 {
397  octave_idx_type n = a.rows ();
398 
399  if (n != a.cols ())
400  {
401  (*current_liboctave_error_handler) ("EIG requires square matrix");
402  return -1;
403  }
404 
405  octave_idx_type info = 0;
406 
407  FloatComplexMatrix atmp = a;
408  FloatComplex *tmp_data = atmp.fortran_vec ();
409 
410  FloatColumnVector wr (n);
411  float *pwr = wr.fortran_vec ();
412 
413  octave_idx_type lwork = -1;
414  FloatComplex dummy_work;
415 
416  octave_idx_type lrwork = 3*n;
417  Array<float> rwork (dim_vector (lrwork, 1));
418  float *prwork = rwork.fortran_vec ();
419 
420  F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
421  F77_CONST_CHAR_ARG2 ("U", 1),
422  n, tmp_data, n, pwr, &dummy_work, lwork,
423  prwork, info
424  F77_CHAR_ARG_LEN (1)
425  F77_CHAR_ARG_LEN (1)));
426 
427  if (info == 0)
428  {
429  lwork = static_cast<octave_idx_type> (dummy_work.real ());
430  Array<FloatComplex> work (dim_vector (lwork, 1));
431  FloatComplex *pwork = work.fortran_vec ();
432 
433  F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
434  F77_CONST_CHAR_ARG2 ("U", 1),
435  n, tmp_data, n, pwr, pwork, lwork, prwork, info
436  F77_CHAR_ARG_LEN (1)
437  F77_CHAR_ARG_LEN (1)));
438 
439  if (info < 0)
440  {
441  (*current_liboctave_error_handler) ("unrecoverable error in cheev");
442  return info;
443  }
444 
445  if (info > 0)
446  {
447  (*current_liboctave_error_handler) ("cheev failed to converge");
448  return info;
449  }
450 
452  v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
453  }
454  else
455  (*current_liboctave_error_handler) ("cheev workspace query failed");
456 
457  return info;
458 }
459 
461 FloatEIG::init (const FloatMatrix& a, const FloatMatrix& b, bool calc_ev)
462 {
464  {
465  (*current_liboctave_error_handler)
466  ("EIG: matrix contains Inf or NaN values");
467  return -1;
468  }
469 
470  octave_idx_type n = a.rows ();
471  octave_idx_type nb = b.rows ();
472 
473  if (n != a.cols () || nb != b.cols ())
474  {
475  (*current_liboctave_error_handler) ("EIG requires square matrix");
476  return -1;
477  }
478 
479  if (n != nb)
480  {
481  (*current_liboctave_error_handler) ("EIG requires same size matrices");
482  return -1;
483  }
484 
485  octave_idx_type info = 0;
486 
487  FloatMatrix tmp = b;
488  float *tmp_data = tmp.fortran_vec ();
489 
490  F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
491  n, tmp_data, n,
492  info
493  F77_CHAR_ARG_LEN (1)
494  F77_CHAR_ARG_LEN (1)));
495 
496  if (a.is_symmetric () && b.is_symmetric () && info == 0)
497  return symmetric_init (a, b, calc_ev);
498 
499  FloatMatrix atmp = a;
500  float *atmp_data = atmp.fortran_vec ();
501 
502  FloatMatrix btmp = b;
503  float *btmp_data = btmp.fortran_vec ();
504 
505  Array<float> ar (dim_vector (n, 1));
506  float *par = ar.fortran_vec ();
507 
508  Array<float> ai (dim_vector (n, 1));
509  float *pai = ai.fortran_vec ();
510 
511  Array<float> beta (dim_vector (n, 1));
512  float *pbeta = beta.fortran_vec ();
513 
514  volatile octave_idx_type nvr = calc_ev ? n : 0;
515  FloatMatrix vr (nvr, nvr);
516  float *pvr = vr.fortran_vec ();
517 
518  octave_idx_type lwork = -1;
519  float dummy_work;
520 
521  float *dummy = 0;
522  octave_idx_type idummy = 1;
523 
524  F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
525  F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
526  n, atmp_data, n, btmp_data, n,
527  par, pai, pbeta,
528  dummy, idummy, pvr, n,
529  &dummy_work, lwork, info
530  F77_CHAR_ARG_LEN (1)
531  F77_CHAR_ARG_LEN (1)));
532 
533  if (info == 0)
534  {
535  lwork = static_cast<octave_idx_type> (dummy_work);
536  Array<float> work (dim_vector (lwork, 1));
537  float *pwork = work.fortran_vec ();
538 
539  F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
540  F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
541  n, atmp_data, n, btmp_data, n,
542  par, pai, pbeta,
543  dummy, idummy, pvr, n,
544  pwork, lwork, info
545  F77_CHAR_ARG_LEN (1)
546  F77_CHAR_ARG_LEN (1)));
547 
548  if (info < 0)
549  {
550  (*current_liboctave_error_handler) ("unrecoverable error in sggev");
551  return info;
552  }
553 
554  if (info > 0)
555  {
556  (*current_liboctave_error_handler) ("sggev failed to converge");
557  return info;
558  }
559 
560  lambda.resize (n);
561  v.resize (nvr, nvr);
562 
563  for (octave_idx_type j = 0; j < n; j++)
564  {
565  if (ai.elem (j) == 0.0)
566  {
567  lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j));
568  for (octave_idx_type i = 0; i < nvr; i++)
569  v.elem (i, j) = vr.elem (i, j);
570  }
571  else
572  {
573  if (j+1 >= n)
574  {
575  (*current_liboctave_error_handler) ("EIG: internal error");
576  return -1;
577  }
578 
579  lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j),
580  ai.elem (j) / beta.elem (j));
581  lambda.elem (j+1) = FloatComplex (ar.elem (j+1) / beta.elem (j+1),
582  ai.elem (j+1) / beta.elem (j+1));
583 
584  for (octave_idx_type i = 0; i < nvr; i++)
585  {
586  float real_part = vr.elem (i, j);
587  float imag_part = vr.elem (i, j+1);
588  v.elem (i, j) = FloatComplex (real_part, imag_part);
589  v.elem (i, j+1) = FloatComplex (real_part, -imag_part);
590  }
591  j++;
592  }
593  }
594  }
595  else
596  (*current_liboctave_error_handler) ("sggev workspace query failed");
597 
598  return info;
599 }
600 
603  bool calc_ev)
604 {
605  octave_idx_type n = a.rows ();
606  octave_idx_type nb = b.rows ();
607 
608  if (n != a.cols () || nb != b.cols ())
609  {
610  (*current_liboctave_error_handler) ("EIG requires square matrix");
611  return -1;
612  }
613 
614  if (n != nb)
615  {
616  (*current_liboctave_error_handler) ("EIG requires same size matrices");
617  return -1;
618  }
619 
620  octave_idx_type info = 0;
621 
622  FloatMatrix atmp = a;
623  float *atmp_data = atmp.fortran_vec ();
624 
625  FloatMatrix btmp = b;
626  float *btmp_data = btmp.fortran_vec ();
627 
628  FloatColumnVector wr (n);
629  float *pwr = wr.fortran_vec ();
630 
631  octave_idx_type lwork = -1;
632  float dummy_work;
633 
634  F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
635  F77_CONST_CHAR_ARG2 ("U", 1),
636  n, atmp_data, n,
637  btmp_data, n,
638  pwr, &dummy_work, lwork, info
639  F77_CHAR_ARG_LEN (1)
640  F77_CHAR_ARG_LEN (1)));
641 
642  if (info == 0)
643  {
644  lwork = static_cast<octave_idx_type> (dummy_work);
645  Array<float> work (dim_vector (lwork, 1));
646  float *pwork = work.fortran_vec ();
647 
648  F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
649  F77_CONST_CHAR_ARG2 ("U", 1),
650  n, atmp_data, n,
651  btmp_data, n,
652  pwr, pwork, lwork, info
653  F77_CHAR_ARG_LEN (1)
654  F77_CHAR_ARG_LEN (1)));
655 
656  if (info < 0)
657  {
658  (*current_liboctave_error_handler) ("unrecoverable error in ssygv");
659  return info;
660  }
661 
662  if (info > 0)
663  {
664  (*current_liboctave_error_handler) ("ssygv failed to converge");
665  return info;
666  }
667 
669  v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
670  }
671  else
672  (*current_liboctave_error_handler) ("ssygv workspace query failed");
673 
674  return info;
675 }
676 
679  bool calc_ev)
680 {
682  {
683  (*current_liboctave_error_handler)
684  ("EIG: matrix contains Inf or NaN values");
685  return -1;
686  }
687 
688  octave_idx_type n = a.rows ();
689  octave_idx_type nb = b.rows ();
690 
691  if (n != a.cols () || nb != b.cols ())
692  {
693  (*current_liboctave_error_handler) ("EIG requires square matrix");
694  return -1;
695  }
696 
697  if (n != nb)
698  {
699  (*current_liboctave_error_handler) ("EIG requires same size matrices");
700  return -1;
701  }
702 
703  octave_idx_type info = 0;
704 
705  FloatComplexMatrix tmp = b;
706  FloatComplex *tmp_data = tmp.fortran_vec ();
707 
708  F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
709  n, tmp_data, n,
710  info
711  F77_CHAR_ARG_LEN (1)
712  F77_CHAR_ARG_LEN (1)));
713 
714  if (a.is_hermitian () && b.is_hermitian () && info == 0)
715  return hermitian_init (a, calc_ev);
716 
717  FloatComplexMatrix atmp = a;
718  FloatComplex *atmp_data = atmp.fortran_vec ();
719 
720  FloatComplexMatrix btmp = b;
721  FloatComplex *btmp_data = btmp.fortran_vec ();
722 
723  FloatComplexColumnVector alpha (n);
724  FloatComplex *palpha = alpha.fortran_vec ();
725 
726  FloatComplexColumnVector beta (n);
727  FloatComplex *pbeta = beta.fortran_vec ();
728 
729  octave_idx_type nvr = calc_ev ? n : 0;
730  FloatComplexMatrix vtmp (nvr, nvr);
731  FloatComplex *pv = vtmp.fortran_vec ();
732 
733  octave_idx_type lwork = -1;
734  FloatComplex dummy_work;
735 
736  octave_idx_type lrwork = 8*n;
737  Array<float> rwork (dim_vector (lrwork, 1));
738  float *prwork = rwork.fortran_vec ();
739 
740  FloatComplex *dummy = 0;
741  octave_idx_type idummy = 1;
742 
743  F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
744  F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
745  n, atmp_data, n, btmp_data, n,
746  palpha, pbeta, dummy, idummy,
747  pv, n, &dummy_work, lwork, prwork, info
748  F77_CHAR_ARG_LEN (1)
749  F77_CHAR_ARG_LEN (1)));
750 
751  if (info == 0)
752  {
753  lwork = static_cast<octave_idx_type> (dummy_work.real ());
754  Array<FloatComplex> work (dim_vector (lwork, 1));
755  FloatComplex *pwork = work.fortran_vec ();
756 
757  F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
758  F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
759  n, atmp_data, n, btmp_data, n,
760  palpha, pbeta, dummy, idummy,
761  pv, n, pwork, lwork, prwork, info
762  F77_CHAR_ARG_LEN (1)
763  F77_CHAR_ARG_LEN (1)));
764 
765  if (info < 0)
766  {
767  (*current_liboctave_error_handler) ("unrecoverable error in cggev");
768  return info;
769  }
770 
771  if (info > 0)
772  {
773  (*current_liboctave_error_handler) ("cggev failed to converge");
774  return info;
775  }
776 
777  lambda.resize (n);
778 
779  for (octave_idx_type j = 0; j < n; j++)
780  lambda.elem (j) = alpha.elem (j) / beta.elem (j);
781 
782  v = vtmp;
783  }
784  else
785  (*current_liboctave_error_handler) ("cggev workspace query failed");
786 
787  return info;
788 }
789 
792  const FloatComplexMatrix& b, bool calc_ev)
793 {
794  octave_idx_type n = a.rows ();
795  octave_idx_type nb = b.rows ();
796 
797  if (n != a.cols () || nb != b.cols ())
798  {
799  (*current_liboctave_error_handler) ("EIG requires square matrix");
800  return -1;
801  }
802 
803  if (n != nb)
804  {
805  (*current_liboctave_error_handler) ("EIG requires same size matrices");
806  return -1;
807  }
808 
809  octave_idx_type info = 0;
810 
811  FloatComplexMatrix atmp = a;
812  FloatComplex *atmp_data = atmp.fortran_vec ();
813 
814  FloatComplexMatrix btmp = b;
815  FloatComplex *btmp_data = btmp.fortran_vec ();
816 
817  FloatColumnVector wr (n);
818  float *pwr = wr.fortran_vec ();
819 
820  octave_idx_type lwork = -1;
821  FloatComplex dummy_work;
822 
823  octave_idx_type lrwork = 3*n;
824  Array<float> rwork (dim_vector (lrwork, 1));
825  float *prwork = rwork.fortran_vec ();
826 
827  F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
828  F77_CONST_CHAR_ARG2 ("U", 1),
829  n, atmp_data, n,
830  btmp_data, n,
831  pwr, &dummy_work, lwork,
832  prwork, info
833  F77_CHAR_ARG_LEN (1)
834  F77_CHAR_ARG_LEN (1)));
835 
836  if (info == 0)
837  {
838  lwork = static_cast<octave_idx_type> (dummy_work.real ());
839  Array<FloatComplex> work (dim_vector (lwork, 1));
840  FloatComplex *pwork = work.fortran_vec ();
841 
842  F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
843  F77_CONST_CHAR_ARG2 ("U", 1),
844  n, atmp_data, n,
845  btmp_data, n,
846  pwr, pwork, lwork, prwork, info
847  F77_CHAR_ARG_LEN (1)
848  F77_CHAR_ARG_LEN (1)));
849 
850  if (info < 0)
851  {
852  (*current_liboctave_error_handler) ("unrecoverable error in zhegv");
853  return info;
854  }
855 
856  if (info > 0)
857  {
858  (*current_liboctave_error_handler) ("zhegv failed to converge");
859  return info;
860  }
861 
863  v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
864  }
865  else
866  (*current_liboctave_error_handler) ("zhegv workspace query failed");
867 
868  return info;
869 }