54 #if defined (DEBUG) || defined (DEBUG_SORT)
63 const double&
BETA,
const double& S,
230 const octave_idx_type&,
const double*,
231 const octave_idx_type&,
double*,
double&
245 const double& beta,
const double& s,
const double&)
248 return (alpha * beta >= 0 ? 1 : -1);
250 return (s >= 0 ? 1 : -1);
255 const double& beta,
const double&,
const double& p)
260 retval = (fabs (alpha) < fabs (beta) ? 1 : -1);
262 retval = (fabs (p) < 1 ? 1 : -1);
265 std::cout <<
"qz: fin: retval=" << retval << std::endl;
273 const double& beta,
const double& s,
const double&)
276 return (alpha * beta < 0 ? 1 : -1);
278 return (s < 0 ? 1 : -1);
283 const double& beta,
const double&,
const double& p)
286 return (fabs (alpha) >= fabs (beta) ? 1 : -1);
288 return (fabs (p) >= 1 ? 1 : -1);
294 DEFUN (qz, args, nargout,
296 @deftypefn {Built-in Function} {@var{lambda} =} qz (@var{A}, @var{B})\n\
297 @deftypefnx {Built-in Function} {@var{lambda} =} qz (@var{A}, @var{B}, @var{opt})\n\
298 QZ@tie{}decomposition of the generalized eigenvalue problem\n\
299 (@math{A x = s B x}).\n\
301 There are three ways to call this function:\n\
303 @item @code{@var{lambda} = qz (@var{A}, @var{B})}\n\
305 Computes the generalized eigenvalues\n\
312 of @math{(A - s B)}.\n\
314 @item @code{[AA, BB, Q, Z, V, W, @var{lambda}] = qz (@var{A}, @var{B})}\n\
316 Computes QZ@tie{}decomposition, generalized eigenvectors, and generalized\n\
317 eigenvalues of @math{(A - s B)}\n\
319 $$ AV = BV{ \\rm diag }(\\lambda) $$\n\
320 $$ W^T A = { \\rm diag }(\\lambda)W^T B $$\n\
321 $$ AA = Q^T AZ, BB = Q^T BZ $$\n\
328 A * V = B * V * diag (@var{lambda})\n\
329 W' * A = diag (@var{lambda}) * W' * B\n\
330 AA = Q * A * Z, BB = Q * B * Z\n\
336 with @var{Q} and @var{Z} orthogonal (unitary)= @var{I}\n\
338 @item @code{[AA,BB,Z@{, @var{lambda}@}] = qz (@var{A}, @var{B}, @var{opt})}\n\
340 As in form [2], but allows ordering of generalized eigenpairs for, e.g.,\n\
341 solution of discrete time algebraic Riccati equations. Form 3 is not\n\
342 available for complex matrices, and does not compute the generalized\n\
343 eigenvectors @var{V}, @var{W}, nor the orthogonal matrix @var{Q}.\n\
347 for ordering eigenvalues of the @nospell{GEP} pencil. The leading block of\n\
348 the revised pencil contains all eigenvalues that satisfy:\n\
351 @item @qcode{\"N\"}\n\
352 = unordered (default)\n\
354 @item @qcode{\"S\"}\n\
355 = small: leading block has all |lambda| @leq{} 1\n\
357 @item @qcode{\"B\"}\n\
358 = big: leading block has all |lambda| @geq{} 1\n\
360 @item @qcode{\"-\"}\n\
361 = negative real part: leading block has all eigenvalues\n\
362 in the open left half-plane\n\
364 @item @qcode{\"+\"}\n\
365 = non-negative real part: leading block has all eigenvalues\n\
366 in the closed right half-plane\n\
371 Note: @code{qz} performs permutation balancing, but not scaling\n\
372 (@pxref{XREFbalance}). The order of output arguments was selected for\n\
373 compatibility with @sc{matlab}.\n\
374 @seealso{eig, balance, lu, chol, hess, qr, qzhess, schur, svd}\n\
378 int nargin = args.
length ();
381 std::cout <<
"qz: nargin = " << nargin
382 <<
", nargout = " << nargout << std::endl;
385 if (nargin < 2 || nargin > 3 || nargout > 7)
390 else if (nargin == 3 && (nargout < 3 || nargout > 4))
392 error (
"qz: invalid number of output arguments for form [3] call");
397 std::cout <<
"qz: determine ordering option" << std::endl;
401 volatile char ord_job = 0;
402 static double safmin;
406 else if (! args(2).is_string ())
408 error (
"qz: OPT must be a string");
413 std::string tmp = args(2).string_value ();
418 if (! (ord_job ==
'N' || ord_job ==
'n'
419 || ord_job ==
'S' || ord_job ==
's'
420 || ord_job ==
'B' || ord_job ==
'b'
421 || ord_job ==
'+' || ord_job ==
'-'))
423 error (
"qz: invalid order option");
433 std::cout <<
"qz: initial value of safmin="
434 << setiosflags (std::ios::scientific)
435 << safmin << std::endl;
443 std::cout <<
"qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl;
451 std::cout <<
"qz: safmin set to "
452 << setiosflags (std::ios::scientific)
453 << safmin << std::endl;
459 std::cout <<
"qz: check argument 1" << std::endl;
466 std::cout <<
"argument 1 dimensions: ("
467 << nn <<
"," << args(0).columns () <<
")"
471 int arg_is_empty =
empty_arg (
"qz", nn, args(0).columns ());
473 if (arg_is_empty < 0)
478 else if (arg_is_empty > 0)
483 else if (args(0).columns () != nn)
493 if (args(0).is_complex_type ())
494 caa = args(0).complex_matrix_value ();
496 aa = args(0).matrix_value ();
502 std::cout <<
"qz: check argument 2" << std::endl;
506 if ((nn != args(1).columns ()) || (nn != args(1).rows ()))
515 if (args(1).is_complex_type ())
516 cbb = args(1).complex_matrix_value ();
518 bb = args(1).matrix_value ();
527 volatile int complex_case
528 = (args(0).is_complex_type () || args(1).is_complex_type ());
530 if (nargin == 3 && complex_case)
532 error (
"qz: cannot re-order complex qz decomposition");
537 Matrix QQ(nn,nn), ZZ(nn,nn),
VR(nn,nn),
VL(nn,nn);
538 RowVector alphar(nn), alphai(nn), betar(nn);
542 char compq = (nargout >= 3 ?
'V' :
'N');
543 char compz = ((nargout >= 4 || nargin == 3)?
'V' :
'N');
546 if (compq ==
'V' || compz ==
'V')
551 QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0);
555 const char bal_job =
'P';
556 RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn);
562 std::cout <<
"qz: performing balancing; CQ=" << std::endl
565 if (args(0).is_real_type ())
568 if (args(1).is_real_type ())
580 nn, ilo, ihi, lscale.fortran_vec (),
581 rscale.fortran_vec (), work.fortran_vec (), info
588 std::cout <<
"qz: performing balancing; QQ=" << std::endl
595 nn, ilo, ihi, lscale.fortran_vec (),
596 rscale.fortran_vec (), work.fortran_vec (), info
609 nn, ilo, ihi, lscale.data (), rscale.data (),
610 nn, QQ.fortran_vec (),
nn, info
616 std::cout <<
"qz: balancing done; QQ=" << std::endl << QQ << std::endl;
626 nn, ilo, ihi, lscale.data (), rscale.data (),
627 nn, ZZ.fortran_vec (),
nn, info
633 std::cout <<
"qz: balancing done; ZZ=" << std::endl << ZZ << std::endl;
639 qz_job = (nargout < 2 ?
'E' :
'S');
658 CZ.fortran_vec (),
nn, info
672 CQ.fortran_vec (),
nn,
673 CZ.fortran_vec (),
nn,
685 nn, ilo, ihi, lscale.data (), rscale.data (),
686 nn,
CQ.fortran_vec (),
nn, info
697 nn, ilo, ihi, lscale.data (), rscale.data (),
698 nn,
CZ.fortran_vec (),
nn, info
707 std::cout <<
"qz: peforming qr decomposition of bb" << std::endl;
714 std::cout <<
"qz: qr (bb) done; now peforming qz decomposition"
721 std::cout <<
"qz: extracted bb" << std::endl;
727 std::cout <<
"qz: updated aa " << std::endl;
728 std::cout <<
"bqr.Q () = " << std::endl << bqr.
Q () << std::endl;
731 std::cout <<
"QQ =" << QQ << std::endl;
738 std::cout <<
"qz: precursors done..." << std::endl;
742 std::cout <<
"qz: compq = " << compq <<
", compz = " << compz
752 ZZ.fortran_vec (),
nn, info
765 nn, alphar.fortran_vec (), alphai.fortran_vec (),
767 ZZ.fortran_vec (),
nn, work.fortran_vec (),
nn, info
777 nn, ilo, ihi, lscale.data (), rscale.data (),
778 nn, QQ.fortran_vec (),
nn, info
784 std::cout <<
"qz: balancing done; QQ=" << std::endl
795 nn, ilo, ihi, lscale.data (), rscale.data (),
796 nn, ZZ.fortran_vec (),
nn, info
802 std::cout <<
"qz: balancing done; ZZ=" << std::endl
810 if (! (ord_job ==
'N' || ord_job ==
'n'))
815 error (
"qz: cannot re-order complex qz decomposition");
821 std::cout <<
"qz: ordering eigenvalues: ord_job = "
822 << ord_job << std::endl;
861 nn, nn, aa.
data (),
nn, work.fortran_vec (), inf_norm
864 double eps = std::numeric_limits<double>::epsilon () * inf_norm *
nn;
867 std::cout <<
"qz: calling dsubsp: aa=" << std::endl;
869 std::cout << std::endl <<
"bb=" << std::endl;
873 std::cout << std::endl <<
"ZZ=" << std::endl;
876 std::cout << std::endl;
877 std::cout <<
"alphar = " << std::endl;
879 std::cout << std::endl <<
"alphai = " << std::endl;
881 std::cout << std::endl <<
"beta = " << std::endl;
883 std::cout << std::endl;
890 ZZ.fortran_vec (), sort_test,
eps, ndim, fail,
894 std::cout <<
"qz: back from dsubsp: aa=" << std::endl;
896 std::cout << std::endl <<
"bb=" << std::endl;
900 std::cout << std::endl <<
"ZZ=" << std::endl;
903 std::cout << std::endl;
913 std::cout <<
"computing gen eig #" << jj << std::endl;
921 else if (aa(jj+1,jj) == 0)
929 std::cout <<
" single gen eig:" << std::endl;
930 std::cout <<
" alphar(" << jj <<
") = " << aa(jj,jj)
932 std::cout <<
" betar(" << jj <<
") = " << bb(jj,jj)
934 std::cout <<
" alphai(" << jj <<
") = 0" << std::endl;
937 alphar(jj) = aa(jj,jj);
939 betar(jj) = bb(jj,jj);
945 std::cout <<
"qz: calling dlag2:" << std::endl;
946 std::cout <<
"safmin="
947 << setiosflags (std::ios::scientific)
948 << safmin << std::endl;
950 for (
int idr = jj; idr <= jj+1; idr++)
952 for (
int idc = jj; idc <= jj+1; idc++)
954 std::cout <<
"aa(" << idr <<
"," << idc <<
")="
955 << aa(idr,idc) << std::endl;
956 std::cout <<
"bb(" << idr <<
"," << idc <<
")="
957 << bb(idr,idc) << std::endl;
965 double scale1, scale2, wr1, wr2,
wi;
966 const double *aa_ptr = aa.
data () + jj * nn + jj;
967 const double *bb_ptr = bb.
data () + jj * nn + jj;
969 (aa_ptr, nn, bb_ptr, nn, safmin,
970 scale1, scale2, wr1, wr2, wi));
973 std::cout <<
"dlag2 returns: scale1=" << scale1
974 <<
"\tscale2=" << scale2 << std::endl
975 <<
"\twr1=" << wr1 <<
"\twr2=" << wr2
976 <<
"\twi=" << wi << std::endl;
987 betar(jj+1) = scale2;
991 alphar(jj) = alphar(jj+1) = wr1;
992 alphai(jj) = -(alphai(jj+1) =
wi);
993 betar(jj) = betar(jj+1) = scale1;
1002 std::cout <<
"qz: back from dsubsp: aa=" << std::endl;
1004 std::cout << std::endl <<
"bb=" << std::endl;
1009 std::cout << std::endl <<
"ZZ=" << std::endl;
1012 std::cout << std::endl <<
"qz: ndim=" << ndim << std::endl
1013 <<
"fail=" << fail << std::endl;
1014 std::cout <<
"alphar = " << std::endl;
1016 std::cout << std::endl <<
"alphai = " << std::endl;
1018 std::cout << std::endl <<
"beta = " << std::endl;
1020 std::cout << std::endl;
1028 if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4))
1034 for (
int ii = 0; ii <
nn; ii++)
1040 for (
int ii = 0; ii <
nn; ii++)
1041 tmp(cnt++) = xalpha(ii) / xbeta(ii);
1048 std::cout <<
"qz: computing generalized eigenvalues" << std::endl;
1054 for (
int ii = 0; ii <
nn; ii++)
1061 for (
int ii = 0; ii <
nn; ii++)
1063 tmp(cnt++) =
Complex(alphar(ii), alphai(ii))/betar(ii);
1073 char side = (nargout == 5 ?
'R' :
'B');
1099 std::cout <<
"qz: computing generalized eigenvectors" << std::endl;
1111 m, work.fortran_vec (), info
1130 else if (aa(jj+1,jj) == 0)
1136 for (
int ii = 0; ii <
nn; ii++)
1137 CVR(ii,jj) =
VR(ii,jj);
1140 for (
int ii = 0; ii <
nn; ii++)
1141 CVL(ii,jj) =
VL(ii,jj);
1147 for (
int ii = 0; ii <
nn; ii++)
1154 for (
int ii = 0; ii <
nn; ii++)
1184 std::cout <<
"qz: sort: retval(3) = gev = " << std::endl;
1186 std::cout << std::endl;
1209 retval(2) =
CQ.hermitian ();
1211 retval(2) = QQ.transpose ();
1219 std::cout <<
"qz: retval(1) = cbb = " << std::endl;
1221 std::cout << std::endl <<
"qz: retval(0) = caa = " <<std::endl;
1223 std::cout << std::endl;
1231 std::cout <<
"qz: retval(1) = bb = " << std::endl;
1233 std::cout << std::endl <<
"qz: retval(0) = aa = " <<std::endl;
1235 std::cout << std::endl;
1247 std::cout <<
"qz: retval(0) = gev = " << gev << std::endl;
1253 error (
"qz: too many return arguments");
1258 std::cout <<
"qz: exiting (at long last)" << std::endl;
F77_RET_T const octave_idx_type const double const octave_idx_type const double double & SCALE1
#define F77_CHAR_ARG_LEN(l)
subroutine xdlamch(cmach, retval)
F77_RET_T octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double const octave_idx_type double * VR
void gripe_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type Complex const octave_idx_type Complex const octave_idx_type Complex Complex Complex const octave_idx_type Complex const octave_idx_type Complex const octave_idx_type double * RWORK
F77_RET_T octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double const octave_idx_type double const octave_idx_type & LDVR
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type Complex const octave_idx_type Complex const octave_idx_type Complex Complex Complex * CQ
F77_RET_T const octave_idx_type double double double const double & EPS
F77_RET_T octave_idx_type * SELECT
F77_RET_T octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double const octave_idx_type double const octave_idx_type const octave_idx_type & MM
OCTINTERP_API void print_usage(void)
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double double double double const octave_idx_type double const octave_idx_type double const octave_idx_type & LWORK
static octave_idx_type nn
octave_idx_type length(void) const
subroutine xdlange(norm, m, n, a, lda, work, retval)
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)
F77_RET_T octave_idx_type const octave_idx_type const Complex const octave_idx_type const Complex const octave_idx_type Complex const octave_idx_type Complex * xVR
#define DEFUN(name, args_name, nargout_name, doc)
void error(const char *fmt,...)
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double double * ALPHAI
int empty_arg(const char *, octave_idx_type nr, octave_idx_type nc)
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double const octave_idx_type double * Z
F77_RET_T F77_FUNC(dggbal, DGGBAL)(F77_CONST_CHAR_ARG_DECL
void gripe_square_matrix_required(const char *name)
#define F77_XFCN(f, F, args)
F77_RET_T const octave_idx_type const double const octave_idx_type const double & SAFMIN
#define F77_CONST_CHAR_ARG2(x, l)
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double const octave_idx_type double const octave_idx_type & LDZ
F77_RET_T const octave_idx_type double const octave_idx_type double const octave_idx_type octave_idx_type octave_idx_type double double double octave_idx_type &INFO F77_CHAR_ARG_LEN_DECL
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type Complex const octave_idx_type Complex const octave_idx_type Complex * ALPHA
F77_RET_T const octave_idx_type double const octave_idx_type double const octave_idx_type octave_idx_type & ILO
ComplexMatrix hermitian(void) const
F77_RET_T const octave_idx_type double const octave_idx_type double const octave_idx_type & LDB
octave_idx_type(* sort_function)(const octave_idx_type &LSIZE, const double &ALPHA, const double &BETA, const double &S, const double &P)
static octave_idx_type fcrhp(const octave_idx_type &lsize, const double &alpha, const double &beta, const double &s, const double &)
F77_RET_T octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double * VL
const T * data(void) const
F77_RET_T const octave_idx_type const double const octave_idx_type const double double double & SCALE2
F77_RET_T octave_idx_type const octave_idx_type const Complex const octave_idx_type const Complex const octave_idx_type Complex const octave_idx_type Complex const octave_idx_type const octave_idx_type octave_idx_type Complex * CWORK
#define panic_impossible()
subroutine dsubsp(NMAX, N, A, B, Z, FTEST, EPS, NDIM, FAIL, IND)
F77_RET_T const octave_idx_type double const octave_idx_type double const octave_idx_type octave_idx_type octave_idx_type double double * RSCALE
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double double double * BETA
void gripe_empty_arg(const char *name, bool is_error)
F77_RET_T const octave_idx_type const double const octave_idx_type const double double double double double & WR2
F77_RET_T const octave_idx_type double const octave_idx_type double const octave_idx_type octave_idx_type octave_idx_type double double double * WORK
F77_RET_T const octave_idx_type double double double const double octave_idx_type octave_idx_type & FAIL
F77_RET_T const octave_idx_type double double double const double octave_idx_type octave_idx_type octave_idx_type * IND
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type const double const double octave_idx_type double const octave_idx_type & LDV
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type Complex const octave_idx_type Complex const octave_idx_type Complex Complex Complex const octave_idx_type Complex * CZ
static octave_idx_type fout(const octave_idx_type &lsize, const double &alpha, const double &beta, const double &, const double &p)
F77_RET_T octave_idx_type const octave_idx_type const Complex const octave_idx_type const Complex const octave_idx_type Complex * xVL
F77_RET_T const octave_idx_type double const octave_idx_type double * B
F77_RET_T const octave_idx_type double double double const double octave_idx_type & NDIM
F77_RET_T const octave_idx_type double const octave_idx_type double const octave_idx_type octave_idx_type octave_idx_type & IHI
static octave_idx_type fin(const octave_idx_type &lsize, const double &alpha, const double &beta, const double &, const double &p)
F77_RET_T const octave_idx_type double const octave_idx_type double const octave_idx_type octave_idx_type octave_idx_type double * LSCALE
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type const double const double octave_idx_type & M
F77_RET_T const octave_idx_type double * A
F77_RET_T const octave_idx_type & N
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 const octave_idx_type const double const octave_idx_type const double double double double double double & WI
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double * Q
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double * ALPHAR
F77_RET_T const octave_idx_type const double const octave_idx_type const double double double double & WR1
std::complex< double > Complex
const T * fortran_vec(void) const
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type const double const double octave_idx_type double * V
F77_RET_T octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double const octave_idx_type & LDVL
F77_RET_T const octave_idx_type const octave_idx_type const octave_idx_type double const octave_idx_type double const octave_idx_type double const octave_idx_type & LDQ
F77_RET_T const octave_idx_type double const octave_idx_type & LDA
F77_RET_T F77_CONST_CHAR_ARG_DECL