24 #if defined (HAVE_CONFIG_H) 38 if (
a.any_element_is_inf_or_nan ())
40 (
"EIG: matrix contains Inf or NaN values");
45 F77_INT n = octave::to_f77_int (
a.rows ());
49 (*current_liboctave_error_handler) (
"EIG requires square matrix");
60 double *pwi =
wi.fortran_vec ();
62 F77_INT tnvr = (calc_rev ? n : 0);
66 F77_INT tnvl = (calc_lev ? n : 0);
77 double *pscale =
scale.fortran_vec ();
89 F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ?
"B" :
"N", 1),
90 F77_CONST_CHAR_ARG2 (calc_lev ?
"V" :
"N", 1),
91 F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
92 F77_CONST_CHAR_ARG2 (
"N", 1),
93 n, tmp_data, n, pwr, pwi, pvl,
94 n, pvr, n, ilo, ihi, pscale,
95 abnrm, prconde, prcondv, &dummy_work,
96 lwork, &dummy_iwork, info
100 F77_CHAR_ARG_LEN (1)));
103 (*current_liboctave_error_handler) (
"dgeevx workspace query failed");
105 lwork =
static_cast<F77_INT> (dummy_work);
109 F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ?
"B" :
"N", 1),
110 F77_CONST_CHAR_ARG2 (calc_lev ?
"V" :
"N", 1),
111 F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
112 F77_CONST_CHAR_ARG2 (
"N", 1),
113 n, tmp_data, n, pwr, pwi, pvl,
114 n, pvr, n, ilo, ihi, pscale,
115 abnrm, prconde, prcondv, pwork,
116 lwork, &dummy_iwork, info
120 F77_CHAR_ARG_LEN (1)));
123 (*current_liboctave_error_handler) (
"unrecoverable error in dgeevx");
126 (*current_liboctave_error_handler) (
"dgeevx failed to converge");
129 F77_INT nvr = (calc_rev ? n : 0);
131 F77_INT nvl = (calc_lev ? n : 0);
134 for (
F77_INT j = 0; j < n; j++)
136 if (
wi.elem (j) == 0.0)
148 (*current_liboctave_error_handler) (
"EIG: internal error");
155 double real_part = vr.
elem (
i, j);
156 double imag_part = vr.
elem (
i, j+1);
163 double real_part = vl.
elem (
i, j);
164 double imag_part = vl.
elem (
i, j+1);
178 F77_INT n = octave::to_f77_int (
a.rows ());
182 (*current_liboctave_error_handler) (
"EIG requires square matrix");
195 F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
196 F77_CONST_CHAR_ARG2 (
"U", 1),
197 n, tmp_data, n, pwr, &dummy_work, lwork, info
199 F77_CHAR_ARG_LEN (1)));
202 (*current_liboctave_error_handler) (
"dsyev workspace query failed");
204 lwork =
static_cast<F77_INT> (dummy_work);
208 F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
209 F77_CONST_CHAR_ARG2 (
"U", 1),
210 n, tmp_data, n, pwr, pwork, lwork, info
212 F77_CHAR_ARG_LEN (1)));
215 (*current_liboctave_error_handler) (
"unrecoverable error in dsyev");
218 (*current_liboctave_error_handler) (
"dsyev failed to converge");
230 if (
a.any_element_is_inf_or_nan ())
232 (
"EIG: matrix contains Inf or NaN values");
234 if (
a.ishermitian ())
237 F77_INT n = octave::to_f77_int (
a.rows ());
241 (*current_liboctave_error_handler) (
"EIG requires square matrix");
251 F77_INT nvr = (calc_rev ? n : 0);
255 F77_INT nvl = (calc_lev ? n : 0);
270 double *pscale =
scale.fortran_vec ();
280 F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ?
"B" :
"N", 1),
281 F77_CONST_CHAR_ARG2 (calc_lev ?
"V" :
"N", 1),
282 F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
283 F77_CONST_CHAR_ARG2 (
"N", 1),
287 pscale, abnrm, prconde, prcondv,
293 F77_CHAR_ARG_LEN (1)));
296 (*current_liboctave_error_handler) (
"zgeevx workspace query failed");
298 lwork =
static_cast<F77_INT> (dummy_work.real ());
302 F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ?
"B" :
"N", 1),
303 F77_CONST_CHAR_ARG2 (calc_lev ?
"V" :
"N", 1),
304 F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
305 F77_CONST_CHAR_ARG2 (
"N", 1),
309 pscale, abnrm, prconde, prcondv,
314 F77_CHAR_ARG_LEN (1)));
317 (*current_liboctave_error_handler) (
"unrecoverable error in zgeevx");
320 (*current_liboctave_error_handler) (
"zgeevx failed to converge");
332 F77_INT n = octave::to_f77_int (
a.rows ());
336 (*current_liboctave_error_handler) (
"EIG requires square matrix");
353 F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
354 F77_CONST_CHAR_ARG2 (
"U", 1),
359 F77_CHAR_ARG_LEN (1)));
362 (*current_liboctave_error_handler) (
"zheev workspace query failed");
364 lwork =
static_cast<F77_INT> (dummy_work.real ());
368 F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
369 F77_CONST_CHAR_ARG2 (
"U", 1),
373 F77_CHAR_ARG_LEN (1)));
376 (*current_liboctave_error_handler) (
"unrecoverable error in zheev");
379 (*current_liboctave_error_handler) (
"zheev failed to converge");
392 if (
a.any_element_is_inf_or_nan () ||
b.any_element_is_inf_or_nan ())
394 (
"EIG: matrix contains Inf or NaN values");
396 F77_INT n = octave::to_f77_int (
a.rows ());
397 F77_INT nb = octave::to_f77_int (
b.rows ());
403 (*current_liboctave_error_handler) (
"EIG requires square matrix");
406 (*current_liboctave_error_handler) (
"EIG requires same size matrices");
411 double *tmp_data =
tmp.fortran_vec ();
415 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (
"L", 1),
418 F77_CHAR_ARG_LEN (1)));
420 if (
a.issymmetric () &&
b.issymmetric () && info == 0)
431 double *par =
ar.fortran_vec ();
439 F77_INT tnvr = (calc_rev ? n : 0);
443 F77_INT tnvl = (calc_lev ? n : 0);
450 F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ?
"V" :
"N", 1),
451 F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
452 n, atmp_data, n, btmp_data, n,
455 &dummy_work, lwork, info
457 F77_CHAR_ARG_LEN (1)));
460 (*current_liboctave_error_handler) (
"dggev workspace query failed");
462 lwork =
static_cast<F77_INT> (dummy_work);
466 F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ?
"V" :
"N", 1),
467 F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
468 n, atmp_data, n, btmp_data, n,
473 F77_CHAR_ARG_LEN (1)));
476 (*current_liboctave_error_handler) (
"unrecoverable error in dggev");
479 (*current_liboctave_error_handler) (
"dggev failed to converge");
482 F77_INT nvr = (calc_rev ? n : 0);
485 F77_INT nvl = (calc_lev ? n : 0);
488 for (
F77_INT j = 0; j < n; j++)
490 if (ai.
elem (j) == 0.0)
501 (*current_liboctave_error_handler) (
"EIG: internal error");
510 double real_part = vr.
elem (
i, j);
511 double imag_part = vr.
elem (
i, j+1);
517 double real_part = vl.
elem (
i, j);
518 double imag_part = vl.
elem (
i, j+1);
533 F77_INT n = octave::to_f77_int (
a.rows ());
534 F77_INT nb = octave::to_f77_int (
b.rows ());
540 (*current_liboctave_error_handler) (
"EIG requires square matrix");
543 (*current_liboctave_error_handler) (
"EIG requires same size matrices");
559 F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
560 F77_CONST_CHAR_ARG2 (
"U", 1),
563 pwr, &dummy_work, lwork, info
565 F77_CHAR_ARG_LEN (1)));
568 (*current_liboctave_error_handler) (
"dsygv workspace query failed");
570 lwork =
static_cast<F77_INT> (dummy_work);
574 F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
575 F77_CONST_CHAR_ARG2 (
"U", 1),
578 pwr, pwork, lwork, info
580 F77_CHAR_ARG_LEN (1)));
583 (*current_liboctave_error_handler) (
"unrecoverable error in dsygv");
586 (*current_liboctave_error_handler) (
"dsygv failed to converge");
597 bool calc_lev,
bool force_qz)
599 if (
a.any_element_is_inf_or_nan () ||
b.any_element_is_inf_or_nan ())
601 (
"EIG: matrix contains Inf or NaN values");
603 F77_INT n = octave::to_f77_int (
a.rows ());
604 F77_INT nb = octave::to_f77_int (
b.rows ());
610 (*current_liboctave_error_handler) (
"EIG requires square matrix");
613 (*current_liboctave_error_handler) (
"EIG requires same size matrices");
622 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (
"L", 1),
625 F77_CHAR_ARG_LEN (1)));
627 if (
a.ishermitian () &&
b.ishermitian () && info == 0)
643 F77_INT nvr = (calc_rev ? n : 0);
647 F77_INT nvl = (calc_lev ? n : 0);
658 F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 (
"N", 1),
659 F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
669 F77_CHAR_ARG_LEN (1)));
672 (*current_liboctave_error_handler) (
"zggev workspace query failed");
674 lwork =
static_cast<F77_INT> (dummy_work.real ());
678 F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 (
"N", 1),
679 F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
688 F77_CHAR_ARG_LEN (1)));
691 (*current_liboctave_error_handler) (
"unrecoverable error in zggev");
694 (*current_liboctave_error_handler) (
"zggev failed to converge");
698 for (
F77_INT j = 0; j < n; j++)
709 bool calc_rev,
bool calc_lev)
711 F77_INT n = octave::to_f77_int (
a.rows ());
712 F77_INT nb = octave::to_f77_int (
b.rows ());
718 (*current_liboctave_error_handler) (
"EIG requires square matrix");
721 (*current_liboctave_error_handler) (
"EIG requires same size matrices");
741 F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
742 F77_CONST_CHAR_ARG2 (
"U", 1),
748 F77_CHAR_ARG_LEN (1)));
751 (*current_liboctave_error_handler) (
"zhegv workspace query failed");
753 lwork =
static_cast<F77_INT> (dummy_work.real ());
757 F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
758 F77_CONST_CHAR_ARG2 (
"U", 1),
763 F77_CHAR_ARG_LEN (1)));
766 (*current_liboctave_error_handler) (
"unrecoverable error in zhegv");
769 (*current_liboctave_error_handler) (
"zhegv failed to converge");
void resize(octave_idx_type n, const Complex &rfv=Complex(0))
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
#define F77_DBLE_CMPLX_ARG(x)
const T * fortran_vec(void) const
T & elem(octave_idx_type n)
#define F77_XFCN(f, F, args)
octave_idx_type symmetric_init(const Matrix &a, bool calc_rev, bool calc_lev)
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
octave_idx_type hermitian_init(const ComplexMatrix &a, bool calc_rev, bool calc_lev)
octave_idx_type init(const Matrix &a, bool calc_rev, bool calc_lev, bool balance)
ComplexColumnVector lambda
friend class ComplexMatrix
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
octave_f77_int_type F77_INT
void scale(Matrix &m, double x, double y, double z)
std::complex< double > Complex
Vector representing the dimensions (size) of an Array.