32 #if defined (HAVE_CONFIG_H) 39 #if defined (DEBUG) || defined (DEBUG_SORT) || defined (DEBUG_EIG) 41 # if defined (DEBUG_EIG) 55 #if defined (DEBUG) || defined (DEBUG_SORT) 170 int nargin = args.length ();
173 std::cout <<
"qz: nargin = " <<
nargin 174 <<
", nargout = " <<
nargout << std::endl;
180 if (
nargin == 3 && (nargout < 3 || nargout > 4))
181 error (
"qz: invalid number of output arguments for form 3 call");
184 std::cout <<
"qz: determine ordering option" << std::endl;
194 std::string opt = args(2).xstring_value (
"qz: OPT must be a string");
197 error (
"qz: OPT must be a non-empty string");
199 ord_job = std::toupper (opt[0]);
203 if (valid_opts.find_first_of (ord_job) == std::string::npos)
204 error (
"qz: invalid order option '%c'", ord_job);
209 F77_CHAR_ARG_LEN (1));
211 #if defined (DEBUG_EIG) 212 std::cout <<
"qz: initial value of safmin=" 213 << setiosflags (std::ios::scientific)
214 << safmin << std::endl;
221 #if defined (DEBUG_EIG) 222 std::cout <<
"qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl;
227 F77_CHAR_ARG_LEN (1));
229 #if defined (DEBUG_EIG) 230 std::cout <<
"qz: safmin set to " 231 << setiosflags (std::ios::scientific)
232 << safmin << std::endl;
238 std::cout <<
"qz: check matrix A" << std::endl;
246 std::cout <<
"Matrix A dimensions: (" <<
nn <<
',' << nc <<
')' << std::endl;
261 if (args(0).iscomplex ())
262 caa = args(0).complex_matrix_value ();
264 aa = args(0).matrix_value ();
267 std::cout <<
"qz: check matrix B" << std::endl;
280 if (args(1).iscomplex ())
281 cbb = args(1).complex_matrix_value ();
283 bb = args(1).matrix_value ();
287 bool complex_case = (args(0).iscomplex () || args(1).iscomplex ());
289 if (
nargin == 3 && complex_case)
290 error (
"qz: cannot re-order complex qz decomposition");
300 char comp_q = (
nargout >= 3 ?
'V' :
'N');
304 if (comp_q ==
'V' || comp_z ==
'V')
306 double *QQptr = QQ.fortran_vec ();
307 double *ZZptr = ZZ.fortran_vec ();
308 std::fill_n (QQptr, QQ.numel (), 0.0);
309 std::fill_n (ZZptr, ZZ.numel (), 0.0);
318 const char bal_job =
'P';
325 std::cout <<
"qz: performing balancing; CQ=" << std::endl
328 if (args(0).isreal ())
331 if (args(1).isreal ())
341 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
344 nn, ilo, ihi, lscale.fortran_vec (),
345 rscale.fortran_vec (), work.fortran_vec (), info
346 F77_CHAR_ARG_LEN (1)));
352 std::cout <<
"qz: performing balancing; QQ=" << std::endl
357 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
359 nn, ilo, ihi, lscale.fortran_vec (),
360 rscale.fortran_vec (), work.fortran_vec (), info
361 F77_CHAR_ARG_LEN (1)));
373 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
374 F77_CONST_CHAR_ARG2 (
"L", 1),
375 nn, ilo, ihi, lscale.data (), rscale.data (),
376 nn, QQ.fortran_vec (),
nn, info
378 F77_CHAR_ARG_LEN (1)));
382 std::cout <<
"qz: balancing done; QQ=" << std::endl << QQ << std::endl;
390 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
391 F77_CONST_CHAR_ARG2 (
"R", 1),
392 nn, ilo, ihi, lscale.data (), rscale.data (),
393 nn, ZZ.fortran_vec (),
nn, info
395 F77_CHAR_ARG_LEN (1)));
399 std::cout <<
"qz: balancing done; ZZ=" << std::endl << ZZ << std::endl;
404 char qz_job = (
nargout < 2 ?
'E' :
'S');
419 (F77_CONST_CHAR_ARG2 (&comp_q, 1),
420 F77_CONST_CHAR_ARG2 (&comp_z, 1),
426 F77_CHAR_ARG_LEN (1)));
431 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
432 F77_CONST_CHAR_ARG2 (&comp_q, 1),
433 F77_CONST_CHAR_ARG2 (&comp_z, 1),
445 F77_CHAR_ARG_LEN (1)));
451 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
452 F77_CONST_CHAR_ARG2 (
"L", 1),
453 nn, ilo, ihi, lscale.data (), rscale.data (),
456 F77_CHAR_ARG_LEN (1)));
463 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
464 F77_CONST_CHAR_ARG2 (
"R", 1),
465 nn, ilo, ihi, lscale.data (), rscale.data (),
468 F77_CHAR_ARG_LEN (1)));
475 std::cout <<
"qz: performing qr decomposition of bb" << std::endl;
482 std::cout <<
"qz: qr (bb) done; now performing qz decomposition" 489 std::cout <<
"qz: extracted bb" << std::endl;
495 std::cout <<
"qz: updated aa " << std::endl;
496 std::cout <<
"bqr.Q () = " << std::endl << bqr.
Q () << std::endl;
499 std::cout <<
"QQ =" << QQ << std::endl;
506 std::cout <<
"qz: precursors done..." << std::endl;
510 std::cout <<
"qz: comp_q = " << comp_q <<
", comp_z = " << comp_z
516 (F77_CONST_CHAR_ARG2 (&comp_q, 1),
517 F77_CONST_CHAR_ARG2 (&comp_z, 1),
520 ZZ.fortran_vec (),
nn, info
522 F77_CHAR_ARG_LEN (1)));
529 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
530 F77_CONST_CHAR_ARG2 (&comp_q, 1),
531 F77_CONST_CHAR_ARG2 (&comp_z, 1),
533 nn, alphar.fortran_vec (), alphai.fortran_vec (),
535 ZZ.fortran_vec (),
nn, work.fortran_vec (),
nn, info
538 F77_CHAR_ARG_LEN (1)));
543 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
544 F77_CONST_CHAR_ARG2 (
"L", 1),
545 nn, ilo, ihi, lscale.data (), rscale.data (),
546 nn, QQ.fortran_vec (),
nn, info
548 F77_CHAR_ARG_LEN (1)));
552 std::cout <<
"qz: balancing done; QQ=" << std::endl
561 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
562 F77_CONST_CHAR_ARG2 (
"R", 1),
563 nn, ilo, ihi, lscale.data (), rscale.data (),
564 nn, ZZ.fortran_vec (),
nn, info
566 F77_CHAR_ARG_LEN (1)));
570 std::cout <<
"qz: balancing done; ZZ=" << std::endl
582 error (
"qz: cannot re-order complex QZ decomposition");
584 #if defined (DEBUG_SORT) 585 std::cout <<
"qz: ordering eigenvalues: ord_job = " 586 << ord_job << std::endl;
591 for (
int j = 0; j <
nn; j++)
596 select(j) = alphar(j)*alphar(j) + alphai(j)*alphai(j) < betar(j)*betar(j);
600 select(j) = alphar(j)*alphar(j) + alphai(j)*alphai(j) >= betar(j)*betar(j);
604 select(j) = alphar(j) * betar(j) >= 0;
608 select(j) = alphar(j) * betar(j) < 0;
620 F77_INT ijob = 0, mm, lrwork3 = 4*
nn+16, liwork =
nn;
630 alphar.fortran_vec (),
631 alphai.fortran_vec (),
634 ZZ.fortran_vec (),
nn,
642 #if defined (DEBUG_SORT) 643 std::cout <<
"qz: back from dtgsen: aa=" << std::endl;
645 std::cout << std::endl <<
"bb=" << std::endl;
650 std::cout << std::endl <<
"ZZ=" << std::endl;
653 std::cout << std::endl <<
"qz: info=" << info << std::endl;
654 std::cout <<
"alphar = " << std::endl;
656 std::cout << std::endl <<
"alphai = " << std::endl;
658 std::cout << std::endl <<
"beta = " << std::endl;
660 std::cout << std::endl;
674 tmp(
i) = xalpha(
i) / xbeta(
i);
681 std::cout <<
"qz: computing generalized eigenvalues" << std::endl;
702 char side = (
nargout == 5 ?
'R' :
'B');
717 (F77_CONST_CHAR_ARG2 (&side, 1),
718 F77_CONST_CHAR_ARG2 (&howmany, 1),
726 F77_CHAR_ARG_LEN (1)));
731 std::cout <<
"qz: computing generalized eigenvectors" << std::endl;
739 (F77_CONST_CHAR_ARG2 (&side, 1),
740 F77_CONST_CHAR_ARG2 (&howmany, 1),
743 m, work.fortran_vec (), info
745 F77_CHAR_ARG_LEN (1)));
762 else if (aa(j+1,j) == 0)
821 std::cout <<
"qz: sort: retval(3) = gev = " << std::endl;
823 std::cout << std::endl;
849 retval(2) = QQ.transpose ();
858 std::cout <<
"qz: retval(1) = cbb = " << std::endl;
860 std::cout << std::endl <<
"qz: retval(0) = caa = " <<std::endl;
862 std::cout << std::endl;
870 std::cout <<
"qz: retval(1) = bb = " << std::endl;
872 std::cout << std::endl <<
"qz: retval(0) = aa = " <<std::endl;
874 std::cout << std::endl;
885 std::cout <<
"qz: retval(0) = gev = " << gev << std::endl;
891 error (
"qz: too many return arguments");
896 std::cout <<
"qz: exiting (at long last)" << std::endl;
void warn_empty_arg(const char *name)
octave_idx_type rows(void) const
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX F77_DBLE_CMPLX F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX * CZ
subroutine xdlamch(cmach, retval)
OCTINTERP_API void print_usage(void)
#define F77_DBLE_CMPLX_ARG(x)
const T * fortran_vec(void) const
static void transpose(octave_idx_type N, const octave_idx_type *ridx, const octave_idx_type *cidx, octave_idx_type *ridx2, octave_idx_type *cidx2)
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
void error(const char *fmt,...)
octave_idx_type columns(void) const
void err_square_matrix_required(const char *fcn, const char *name)
#define F77_XFCN(f, F, args)
octave_f77_int_type F77_LOGICAL
void octave_print_internal(std::ostream &os, const float_display_format &fmt, bool d, bool pr_as_read_syntax)
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function xample nargout(@histc)
F77_RET_T F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * VL
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
#define panic_impossible()
F77_RET_T F77_FUNC(xerbla, XERBLA)
F77_RET_T F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * VR
N Dimensional Array with copy-on-write semantics.
octave_f77_int_type F77_INT
std::complex< double > Complex
Vector representing the dimensions (size) of an Array.
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX const F77_INT F77_DBLE_CMPLX F77_DBLE_CMPLX F77_DBLE_CMPLX * CQ
If this string is the system will ring the terminal sometimes it is useful to be able to print the original representation of the string
ComplexMatrix hermitian(void) const