31 #if defined (HAVE_CONFIG_H)
53 #if defined (DEBUG) || defined (DEBUG_SORT)
62 const double&
BETA,
const double& S,
87 const double& beta,
const double&
s,
const double&)
90 return (alpha * beta >= 0 ? 1 : -1);
92 return (s >= 0 ? 1 : -1);
97 const double& beta,
const double&,
const double&
p)
102 retval = (fabs (alpha) < fabs (beta) ? 1 : -1);
104 retval = (fabs (p) < 1 ? 1 : -1);
107 std::cout <<
"qz: fin: retval=" << retval << std::endl;
115 const double& beta,
const double&
s,
const double&)
118 return (alpha * beta < 0 ? 1 : -1);
120 return (s < 0 ? 1 : -1);
125 const double& beta,
const double&,
const double&
p)
128 return (fabs (alpha) >= fabs (beta) ? 1 : -1);
130 return (fabs (p) >= 1 ? 1 : -1);
221 std::cout <<
"qz: nargin = " << nargin
222 <<
", nargout = " <<
nargout << std::endl;
225 if (nargin < 2 || nargin > 3 ||
nargout > 7)
228 if (nargin == 3 && (nargout < 3 || nargout > 4))
229 error (
"qz: invalid number of output arguments for form [3] call");
232 std::cout <<
"qz: determine ordering option" << std::endl;
236 volatile char ord_job = 0;
237 static double safmin;
248 if (! (ord_job ==
'N' || ord_job ==
'n'
249 || ord_job ==
'S' || ord_job ==
's'
250 || ord_job ==
'B' || ord_job ==
'b'
251 || ord_job ==
'+' || ord_job ==
'-'))
252 error (
"qz: invalid order option");
257 F77_CHAR_ARG_LEN (1));
259 #if defined (DEBUG_EIG)
260 std::cout <<
"qz: initial value of safmin="
261 << setiosflags (std::ios::scientific)
262 << safmin << std::endl;
269 #if defined (DEBUG_EIG)
270 std::cout <<
"qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl;
275 F77_CHAR_ARG_LEN (1));
277 #if defined (DEBUG_EIG)
278 std::cout <<
"qz: safmin set to "
279 << setiosflags (std::ios::scientific)
280 << safmin << std::endl;
286 std::cout <<
"qz: check argument 1" << std::endl;
293 std::cout <<
"argument 1 dimensions: ("
294 << nn <<
"," <<
args(0).columns () <<
")"
300 if (
args(0).is_empty ())
305 else if (
args(0).columns () != nn)
312 if (
args(0).is_complex_type ())
313 caa =
args(0).complex_matrix_value ();
315 aa =
args(0).matrix_value ();
318 std::cout <<
"qz: check argument 2" << std::endl;
322 if ((nn !=
args(1).columns ()) || (nn !=
args(1).rows ()))
328 if (
args(1).is_complex_type ())
329 cbb =
args(1).complex_matrix_value ();
331 bb =
args(1).matrix_value ();
337 volatile int complex_case
338 = (
args(0).is_complex_type () ||
args(1).is_complex_type ());
340 if (nargin == 3 && complex_case)
341 error (
"qz: cannot re-order complex qz decomposition");
344 Matrix QQ(nn,nn), ZZ(nn,nn),
VR(nn,nn),
VL(nn,nn);
345 RowVector alphar(nn), alphai(nn), betar(nn);
349 char compq = (
nargout >= 3 ?
'V' :
'N');
350 char compz = ((
nargout >= 4 || nargin == 3)?
'V' :
'N');
353 if (compq ==
'V' || compz ==
'V')
358 QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0);
362 const char bal_job =
'P';
363 RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn);
369 std::cout <<
"qz: performing balancing; CQ=" << std::endl
372 if (
args(0).is_real_type ())
375 if (
args(1).is_real_type ())
385 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
388 nn, ilo, ihi, lscale.fortran_vec (),
389 rscale.fortran_vec (), work.fortran_vec (), info
390 F77_CHAR_ARG_LEN (1)));
396 std::cout <<
"qz: performing balancing; QQ=" << std::endl
401 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
403 nn, ilo, ihi, lscale.fortran_vec (),
404 rscale.fortran_vec (), work.fortran_vec (), info
405 F77_CHAR_ARG_LEN (1)));
415 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
416 F77_CONST_CHAR_ARG2 (
"L", 1),
417 nn, ilo, ihi, lscale.data (), rscale.data (),
418 nn, QQ.fortran_vec (),
nn, info
420 F77_CHAR_ARG_LEN (1)));
424 std::cout <<
"qz: balancing done; QQ=" << std::endl << QQ << std::endl;
432 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
433 F77_CONST_CHAR_ARG2 (
"R", 1),
434 nn, ilo, ihi, lscale.data (), rscale.data (),
435 nn, ZZ.fortran_vec (),
nn, info
437 F77_CHAR_ARG_LEN (1)));
441 std::cout <<
"qz: balancing done; ZZ=" << std::endl << ZZ << std::endl;
447 qz_job = (
nargout < 2 ?
'E' :
'S');
462 (F77_CONST_CHAR_ARG2 (&compq, 1),
463 F77_CONST_CHAR_ARG2 (&compz, 1),
469 F77_CHAR_ARG_LEN (1)));
474 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
475 F77_CONST_CHAR_ARG2 (&compq, 1),
476 F77_CONST_CHAR_ARG2 (&compz, 1),
487 F77_CHAR_ARG_LEN (1)));
493 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
494 F77_CONST_CHAR_ARG2 (
"L", 1),
495 nn, ilo, ihi, lscale.data (), rscale.data (),
498 F77_CHAR_ARG_LEN (1)));
505 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
506 F77_CONST_CHAR_ARG2 (
"R", 1),
507 nn, ilo, ihi, lscale.data (), rscale.data (),
510 F77_CHAR_ARG_LEN (1)));
517 std::cout <<
"qz: peforming qr decomposition of bb" << std::endl;
524 std::cout <<
"qz: qr (bb) done; now peforming qz decomposition"
531 std::cout <<
"qz: extracted bb" << std::endl;
537 std::cout <<
"qz: updated aa " << std::endl;
538 std::cout <<
"bqr.Q () = " << std::endl << bqr.
Q () << std::endl;
541 std::cout <<
"QQ =" << QQ << std::endl;
548 std::cout <<
"qz: precursors done..." << std::endl;
552 std::cout <<
"qz: compq = " << compq <<
", compz = " << compz
558 (F77_CONST_CHAR_ARG2 (&compq, 1),
559 F77_CONST_CHAR_ARG2 (&compz, 1),
562 ZZ.fortran_vec (),
nn, info
564 F77_CHAR_ARG_LEN (1)));
571 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
572 F77_CONST_CHAR_ARG2 (&compq, 1),
573 F77_CONST_CHAR_ARG2 (&compz, 1),
575 nn, alphar.fortran_vec (), alphai.fortran_vec (),
577 ZZ.fortran_vec (),
nn, work.fortran_vec (),
nn, info
580 F77_CHAR_ARG_LEN (1)));
585 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
586 F77_CONST_CHAR_ARG2 (
"L", 1),
587 nn, ilo, ihi, lscale.data (), rscale.data (),
588 nn, QQ.fortran_vec (),
nn, info
590 F77_CHAR_ARG_LEN (1)));
594 std::cout <<
"qz: balancing done; QQ=" << std::endl
603 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
604 F77_CONST_CHAR_ARG2 (
"R", 1),
605 nn, ilo, ihi, lscale.data (), rscale.data (),
606 nn, ZZ.fortran_vec (),
nn, info
608 F77_CHAR_ARG_LEN (1)));
612 std::cout <<
"qz: balancing done; ZZ=" << std::endl
620 if (! (ord_job ==
'N' || ord_job ==
'n'))
624 error (
"qz: cannot re-order complex qz decomposition");
626 #if defined (DEBUG_SORT)
627 std::cout <<
"qz: ordering eigenvalues: ord_job = "
628 << ord_job << std::endl;
666 (F77_CONST_CHAR_ARG2 (
"I", 1),
667 nn, nn, aa.
data (),
nn, work.fortran_vec (), inf_norm
668 F77_CHAR_ARG_LEN (1)));
670 double eps = std::numeric_limits<double>::epsilon () * inf_norm *
nn;
672 #if defined (DEBUG_SORT)
673 std::cout <<
"qz: calling dsubsp: aa=" << std::endl;
675 std::cout << std::endl <<
"bb=" << std::endl;
679 std::cout << std::endl <<
"ZZ=" << std::endl;
682 std::cout << std::endl;
683 std::cout <<
"alphar = " << std::endl;
685 std::cout << std::endl <<
"alphai = " << std::endl;
687 std::cout << std::endl <<
"beta = " << std::endl;
689 std::cout << std::endl;
696 ZZ.fortran_vec (), sort_test,
eps, ndim, fail,
700 std::cout <<
"qz: back from dsubsp: aa=" << std::endl;
702 std::cout << std::endl <<
"bb=" << std::endl;
706 std::cout << std::endl <<
"ZZ=" << std::endl;
709 std::cout << std::endl;
718 #if defined (DEBUG_EIG)
719 std::cout <<
"computing gen eig #" << jj << std::endl;
727 else if (aa(jj+1,jj) == 0)
734 #if defined (DEBUG_EIG)
735 std::cout <<
" single gen eig:" << std::endl;
736 std::cout <<
" alphar(" << jj <<
") = " << aa(jj,jj)
738 std::cout <<
" betar(" << jj <<
") = " << bb(jj,jj)
740 std::cout <<
" alphai(" << jj <<
") = 0" << std::endl;
743 alphar(jj) = aa(jj,jj);
745 betar(jj) = bb(jj,jj);
750 #if defined (DEBUG_EIG)
751 std::cout <<
"qz: calling dlag2:" << std::endl;
752 std::cout <<
"safmin="
753 << setiosflags (std::ios::scientific)
754 << safmin << std::endl;
756 for (
int idr = jj; idr <= jj+1; idr++)
758 for (
int idc = jj; idc <= jj+1; idc++)
760 std::cout <<
"aa(" << idr <<
"," << idc <<
")="
761 << aa(idr,idc) << std::endl;
762 std::cout <<
"bb(" << idr <<
"," << idc <<
")="
763 << bb(idr,idc) << std::endl;
771 double scale1, scale2, wr1, wr2,
wi;
772 const double *aa_ptr = aa.
data () + jj * nn + jj;
773 const double *bb_ptr = bb.
data () + jj * nn + jj;
775 (aa_ptr, nn, bb_ptr, nn, safmin,
776 scale1, scale2, wr1, wr2, wi));
778 #if defined (DEBUG_EIG)
779 std::cout <<
"dlag2 returns: scale1=" << scale1
780 <<
"\tscale2=" << scale2 << std::endl
781 <<
"\twr1=" << wr1 <<
"\twr2=" << wr2
782 <<
"\twi=" << wi << std::endl;
793 betar(jj+1) = scale2;
797 alphar(jj) = alphar(jj+1) = wr1;
798 alphai(jj) = -(alphai(jj+1) =
wi);
799 betar(jj) = betar(jj+1) = scale1;
807 #if defined (DEBUG_SORT)
808 std::cout <<
"qz: back from dsubsp: aa=" << std::endl;
810 std::cout << std::endl <<
"bb=" << std::endl;
815 std::cout << std::endl <<
"ZZ=" << std::endl;
818 std::cout << std::endl <<
"qz: ndim=" << ndim << std::endl
819 <<
"fail=" << fail << std::endl;
820 std::cout <<
"alphar = " << std::endl;
822 std::cout << std::endl <<
"alphai = " << std::endl;
824 std::cout << std::endl <<
"beta = " << std::endl;
826 std::cout << std::endl;
839 for (
int ii = 0; ii <
nn; ii++)
845 for (
int ii = 0; ii <
nn; ii++)
846 tmp(cnt++) = xalpha(ii) / xbeta(ii);
853 std::cout <<
"qz: computing generalized eigenvalues" << std::endl;
859 for (
int ii = 0; ii <
nn; ii++)
866 for (
int ii = 0; ii <
nn; ii++)
868 tmp(cnt++) =
Complex(alphar(ii), alphai(ii))/betar(ii);
878 char side = (
nargout == 5 ?
'R' :
'B');
893 (F77_CONST_CHAR_ARG2 (&side, 1),
894 F77_CONST_CHAR_ARG2 (&howmny, 1),
901 F77_CHAR_ARG_LEN (1)));
906 std::cout <<
"qz: computing generalized eigenvectors" << std::endl;
914 (F77_CONST_CHAR_ARG2 (&side, 1),
915 F77_CONST_CHAR_ARG2 (&howmny, 1),
918 m, work.fortran_vec (), info
920 F77_CHAR_ARG_LEN (1)));
937 else if (aa(jj+1,jj) == 0)
943 for (
int ii = 0; ii <
nn; ii++)
944 CVR(ii,jj) =
VR(ii,jj);
947 for (
int ii = 0; ii <
nn; ii++)
948 CVL(ii,jj) =
VL(ii,jj);
954 for (
int ii = 0; ii <
nn; ii++)
961 for (
int ii = 0; ii <
nn; ii++)
991 std::cout <<
"qz: sort: retval(3) = gev = " << std::endl;
993 std::cout << std::endl;
1018 retval(2) = QQ.transpose ();
1026 std::cout <<
"qz: retval(1) = cbb = " << std::endl;
1028 std::cout << std::endl <<
"qz: retval(0) = caa = " <<std::endl;
1030 std::cout << std::endl;
1038 std::cout <<
"qz: retval(1) = bb = " << std::endl;
1040 std::cout << std::endl <<
"qz: retval(0) = aa = " <<std::endl;
1042 std::cout << std::endl;
1053 std::cout <<
"qz: retval(0) = gev = " << gev << std::endl;
1059 error (
"qz: too many return arguments");
1064 std::cout <<
"qz: exiting (at long last)" << std::endl;
void warn_empty_arg(const char *name)
octave_idx_type(* sort_function)(const octave_idx_type &LSIZE, const double &ALPHA, const double &BETA, const double &S, const double &P)
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)
F77_RET_T const F77_INT & N
OCTINTERP_API void print_usage(void)
F77_RET_T const F77_INT F77_DBLE * A
static octave_idx_type nn
F77_RET_T const F77_INT F77_DBLE F77_DBLE F77_DBLE const F77_DBLE F77_INT F77_INT & FAIL
#define F77_DBLE_CMPLX_ARG(x)
subroutine xdlange(norm, m, n, a, lda, work, retval)
F77_RET_T const F77_INT F77_DBLE F77_DBLE * B
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)
void error(const char *fmt,...)
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 * ALPHA
F77_RET_T const F77_INT F77_DBLE F77_DBLE F77_DBLE const F77_DBLE & EPS
void err_square_matrix_required(const char *fcn, const char *name)
#define F77_XFCN(f, F, args)
F77_RET_T const F77_INT F77_DBLE F77_DBLE F77_DBLE * Z
ComplexMatrix hermitian(void) const
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)
nd deftypefn *octave_map m
F77_RET_T F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * VL
static octave_idx_type fcrhp(const octave_idx_type &lsize, const double &alpha, const double &beta, const double &s, const double &)
const T * data(void) const
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
#define panic_impossible()
subroutine dsubsp(NMAX, N, A, B, Z, FTEST, EPS, NDIM, FAIL, IND)
F77_RET_T const F77_INT F77_DBLE F77_DBLE F77_DBLE const F77_DBLE F77_INT F77_INT F77_INT * IND
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE F77_DBLE F77_DBLE * BETA
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
static octave_idx_type fout(const octave_idx_type &lsize, const double &alpha, const double &beta, const double &, const double &p)
static octave_idx_type fin(const octave_idx_type &lsize, const double &alpha, const double &beta, const double &, const double &p)
void octave_print_internal(std::ostream &, char, bool)
static octave_idx_type folhp(const octave_idx_type &lsize, const double &alpha, const double &beta, const double &s, const double &)
F77_RET_T F77_FUNC(dsubsp, DSUBSP)(const F77_INT &NMAX
std::complex< double > Complex
const T * fortran_vec(void) const
F77_RET_T const F77_INT F77_DBLE F77_DBLE F77_DBLE const F77_DBLE F77_INT & NDIM
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