24 #if defined (HAVE_CONFIG_H) 39 if (
a.any_element_is_inf_or_nan ())
41 (
"EIG: matrix contains Inf or NaN values");
46 F77_INT n = octave::to_f77_int (
a.rows ());
50 (*current_liboctave_error_handler) (
"EIG requires square matrix");
61 float *pwi =
wi.fortran_vec ();
63 volatile F77_INT nvr = (calc_rev ? n : 0);
67 volatile F77_INT nvl = (calc_lev ? n : 0);
78 float *pscale =
scale.fortran_vec ();
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,
96 ilo, ihi, pscale, abnrm, prconde, prcondv,
97 &dummy_work, lwork, &dummy_iwork, info
101 F77_CHAR_ARG_LEN (1)));
104 (*current_liboctave_error_handler) (
"sgeevx workspace query failed");
106 lwork =
static_cast<F77_INT> (dummy_work);
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,
116 ilo, ihi, pscale, abnrm, prconde, prcondv,
117 pwork, lwork, &dummy_iwork, info
121 F77_CHAR_ARG_LEN (1)));
124 (*current_liboctave_error_handler) (
"unrecoverable error in sgeevx");
127 (*current_liboctave_error_handler) (
"sgeevx failed to converge");
133 for (
F77_INT j = 0; j < n; j++)
135 if (
wi.elem (j) == 0.0)
147 (*current_liboctave_error_handler) (
"EIG: internal error");
154 float real_part = vr.
elem (
i, j);
155 float imag_part = vr.
elem (
i, j+1);
161 float real_part = vl.
elem (
i, j);
162 float imag_part = vl.
elem (
i, j+1);
176 F77_INT n = octave::to_f77_int (
a.rows ());
180 (*current_liboctave_error_handler) (
"EIG requires square matrix");
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
197 F77_CHAR_ARG_LEN (1)));
200 (*current_liboctave_error_handler) (
"ssyev workspace query failed");
202 lwork =
static_cast<F77_INT> (dummy_work);
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
210 F77_CHAR_ARG_LEN (1)));
213 (*current_liboctave_error_handler) (
"unrecoverable error in ssyev");
216 (*current_liboctave_error_handler) (
"ssyev failed to converge");
229 if (
a.any_element_is_inf_or_nan ())
231 (
"EIG: matrix contains Inf or NaN values");
233 if (
a.ishermitian ())
236 F77_INT n = octave::to_f77_int (
a.rows ());
240 (*current_liboctave_error_handler) (
"EIG requires square matrix");
250 F77_INT nvr = (calc_rev ? n : 0);
254 F77_INT nvl = (calc_lev ? n : 0);
269 float *pscale =
scale.fortran_vec ();
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),
285 ilo, ihi, pscale, abnrm, prconde, prcondv,
290 F77_CHAR_ARG_LEN (1)));
293 (*current_liboctave_error_handler) (
"cgeevx workspace query failed");
295 lwork =
static_cast<F77_INT> (dummy_work.real ());
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),
305 ilo, ihi, pscale, abnrm, prconde, prcondv,
310 F77_CHAR_ARG_LEN (1)));
313 (*current_liboctave_error_handler) (
"unrecoverable error in cgeevx");
316 (*current_liboctave_error_handler) (
"cgeevx failed to converge");
329 F77_INT n = octave::to_f77_int (
a.rows ());
333 (*current_liboctave_error_handler) (
"EIG requires square matrix");
350 F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
351 F77_CONST_CHAR_ARG2 (
"U", 1),
356 F77_CHAR_ARG_LEN (1)));
359 (*current_liboctave_error_handler) (
"cheev workspace query failed");
361 lwork =
static_cast<F77_INT> (dummy_work.real ());
365 F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
366 F77_CONST_CHAR_ARG2 (
"U", 1),
370 F77_CHAR_ARG_LEN (1)));
373 (*current_liboctave_error_handler) (
"unrecoverable error in cheev");
376 (*current_liboctave_error_handler) (
"cheev failed to converge");
387 bool calc_lev,
bool force_qz)
389 if (
a.any_element_is_inf_or_nan () ||
b.any_element_is_inf_or_nan ())
391 (
"EIG: matrix contains Inf or NaN values");
393 F77_INT n = octave::to_f77_int (
a.rows ());
394 F77_INT nb = octave::to_f77_int (
b.rows ());
400 (*current_liboctave_error_handler) (
"EIG requires square matrix");
403 (*current_liboctave_error_handler) (
"EIG requires same size matrices");
408 float *tmp_data =
tmp.fortran_vec ();
411 F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 (
"L", 1),
414 F77_CHAR_ARG_LEN (1)));
416 if (
a.issymmetric () &&
b.issymmetric () && info == 0)
427 float *par =
ar.fortran_vec ();
435 volatile F77_INT nvr = (calc_rev ? n : 0);
439 volatile F77_INT nvl = (calc_lev ? n : 0);
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,
451 &dummy_work, lwork, info
453 F77_CHAR_ARG_LEN (1)));
456 (*current_liboctave_error_handler) (
"sggev workspace query failed");
458 lwork =
static_cast<F77_INT> (dummy_work);
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,
469 F77_CHAR_ARG_LEN (1)));
472 (*current_liboctave_error_handler) (
"unrecoverable error in sggev");
475 (*current_liboctave_error_handler) (
"sggev failed to converge");
482 for (
F77_INT j = 0; j < n; j++)
484 if (ai.
elem (j) == 0.0)
496 (*current_liboctave_error_handler) (
"EIG: internal error");
505 float real_part = vr.
elem (
i, j);
506 float imag_part = vr.
elem (
i, j+1);
512 float real_part = vl.
elem (
i, j);
513 float imag_part = vl.
elem (
i, j+1);
526 bool calc_rev,
bool calc_lev)
528 F77_INT n = octave::to_f77_int (
a.rows ());
529 F77_INT nb = octave::to_f77_int (
b.rows ());
535 (*current_liboctave_error_handler) (
"EIG requires square matrix");
538 (*current_liboctave_error_handler) (
"EIG requires same size matrices");
554 F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
555 F77_CONST_CHAR_ARG2 (
"U", 1),
558 pwr, &dummy_work, lwork, info
560 F77_CHAR_ARG_LEN (1)));
563 (*current_liboctave_error_handler) (
"ssygv workspace query failed");
565 lwork =
static_cast<F77_INT> (dummy_work);
569 F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
570 F77_CONST_CHAR_ARG2 (
"U", 1),
573 pwr, pwork, lwork, info
575 F77_CHAR_ARG_LEN (1)));
578 (*current_liboctave_error_handler) (
"unrecoverable error in ssygv");
581 (*current_liboctave_error_handler) (
"ssygv failed to converge");
592 bool calc_rev,
bool calc_lev,
bool force_qz)
594 if (
a.any_element_is_inf_or_nan () ||
b.any_element_is_inf_or_nan ())
596 (
"EIG: matrix contains Inf or NaN values");
598 F77_INT n = octave::to_f77_int (
a.rows ());
599 F77_INT nb = octave::to_f77_int (
b.rows ());
605 (*current_liboctave_error_handler) (
"EIG requires square matrix");
608 (*current_liboctave_error_handler) (
"EIG requires same size matrices");
617 F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (
"L", 1),
620 F77_CHAR_ARG_LEN (1)));
622 if (
a.ishermitian () &&
b.ishermitian () && info == 0)
638 F77_INT nvr = (calc_rev ? n : 0);
642 F77_INT nvl = (calc_lev ? n : 0);
653 F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ?
"V" :
"N", 1),
654 F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
661 F77_CHAR_ARG_LEN (1)));
664 (*current_liboctave_error_handler) (
"cggev workspace query failed");
666 lwork =
static_cast<F77_INT> (dummy_work.real ());
670 F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ?
"V" :
"N", 1),
671 F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
678 F77_CHAR_ARG_LEN (1)));
681 (*current_liboctave_error_handler) (
"unrecoverable error in cggev");
684 (*current_liboctave_error_handler) (
"cggev failed to converge");
688 for (
F77_INT j = 0; j < n; j++)
700 bool calc_rev,
bool calc_lev)
702 F77_INT n = octave::to_f77_int (
a.rows ());
703 F77_INT nb = octave::to_f77_int (
b.rows ());
709 (*current_liboctave_error_handler) (
"EIG requires square matrix");
712 (*current_liboctave_error_handler) (
"EIG requires same size matrices");
732 F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
733 F77_CONST_CHAR_ARG2 (
"U", 1),
739 F77_CHAR_ARG_LEN (1)));
742 (*current_liboctave_error_handler) (
"zhegv workspace query failed");
744 lwork =
static_cast<F77_INT> (dummy_work.real ());
748 F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ?
"V" :
"N", 1),
749 F77_CONST_CHAR_ARG2 (
"U", 1),
754 F77_CHAR_ARG_LEN (1)));
757 (*current_liboctave_error_handler) (
"unrecoverable error in zhegv");
760 (*current_liboctave_error_handler) (
"zhegv failed to converge");
void resize(octave_idx_type nr, octave_idx_type nc, const FloatComplex &rfv=FloatComplex(0))
void resize(octave_idx_type n, const FloatComplex &rfv=FloatComplex(0))
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
const T * fortran_vec(void) const
T & elem(octave_idx_type n)
#define F77_XFCN(f, F, args)
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 init(const FloatMatrix &a, bool calc_rev, bool calc_lev, bool balance)
FloatComplexColumnVector lambda
friend class FloatComplexMatrix
octave_idx_type hermitian_init(const FloatComplexMatrix &a, bool calc_rev, bool calc_lev)
octave_idx_type symmetric_init(const FloatMatrix &a, bool calc_rev, bool calc_lev)
octave_f77_int_type F77_INT
void scale(Matrix &m, double x, double y, double z)
std::complex< float > FloatComplex
Vector representing the dimensions (size) of an Array.