Mercurial > octave
view libinterp/corefcn/qz.cc @ 33557:bfa4ab12b323 bytecode-interpreter tip
maint: Merge default to bytecode-interpreter
author | Arun Giridhar <arungiridhar@gmail.com> |
---|---|
date | Thu, 09 May 2024 18:18:53 -0400 |
parents | f53ac65ffba6 |
children |
line wrap: on
line source
//////////////////////////////////////////////////////////////////////// // // Copyright (C) 1998-2024 The Octave Project Developers // // See the file COPYRIGHT.md in the top-level directory of this // distribution or <https://octave.org/copyright/>. // // This file is part of Octave. // // Octave is free software: you can redistribute it and/or modify it // under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // Octave is distributed in the hope that it will be useful, but // WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with Octave; see the file COPYING. If not, see // <https://www.gnu.org/licenses/>. // //////////////////////////////////////////////////////////////////////// // Generalized eigenvalue balancing via LAPACK // Originally written by A. S. Hodel <scotte@eng.auburn.edu>, but is // substantially different with the change to use LAPACK. #undef DEBUG #undef DEBUG_SORT #undef DEBUG_EIG #if defined (HAVE_CONFIG_H) # include "config.h" #endif #include <cctype> #include <cmath> #if defined (DEBUG_EIG) # include <iomanip> #endif #include "f77-fcn.h" #include "lo-lapack-proto.h" #include "qr.h" #include "quit.h" #include "defun.h" #include "error.h" #include "errwarn.h" #include "ovl.h" #if defined (DEBUG) || defined (DEBUG_SORT) # include "pager.h" # include "pr-output.h" #endif OCTAVE_BEGIN_NAMESPACE(octave) DEFUN (qz, args, nargout, doc: /* -*- texinfo -*- @deftypefn {} {[@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{V}, @var{W}] =} qz (@var{A}, @var{B}) @deftypefnx {} {[@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{V}, @var{W}] =} qz (@var{A}, @var{B}, @var{opt}) Compute the QZ@tie{}decomposition of a generalized eigenvalue problem. The generalized eigenvalue problem is defined as @tex $$A x = \lambda B x$$ @end tex @ifnottex @math{A x = @var{lambda} B x} @end ifnottex There are two calling forms of the function: @enumerate @item @code{[@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{V}, @var{W}, @var{lambda}] = qz (@var{A}, @var{B})} Compute the complex QZ@tie{}decomposition, generalized eigenvectors, and generalized eigenvalues. @tex $$ AA = Q^T AZ, BB = Q^T BZ $$ $$ { \rm diag }(BB)AV = BV{ \rm diag }(AA) $$ $$ { \rm diag }(BB) W^T A = { \rm diag }(AA) W^T B $$ @end tex @ifnottex @example @group @var{AA} = @var{Q} * @var{A} * @var{Z}, @var{BB} = @var{Q} * @var{B} * @var{Z} @var{A} * @var{V} * diag (diag (@var{BB})) = @var{B} * @var{V} * diag (diag (@var{AA})) diag (diag (@var{BB})) * @var{W}' * @var{A} = diag (diag (@var{AA})) * @var{W}' * @var{B} @end group @end example @end ifnottex with @var{AA} and @var{BB} upper triangular, and @var{Q} and @var{Z} unitary. The matrices @var{V} and @var{W} respectively contain the right and left generalized eigenvectors. @item @code{[@var{AA}, @var{BB}, @var{Z} @{, @var{lambda}@}] = qz (@var{A}, @var{B}, @var{opt})} The @var{opt} argument must be equal to either @qcode{"real"} or @qcode{"complex"}. If it is equal to @qcode{"complex"}, then this calling form is equivalent to the first one with only two input arguments. If @var{opt} is equal to @qcode{"real"}, then the real QZ decomposition is computed. In particular, @var{AA} is only guaranteed to be quasi-upper triangular with 1-by-1 and 2-by-2 blocks on the diagonal, and @var{Q} and @var{Z} are orthogonal. The identities mentioned above for right and left generalized eigenvectors are only verified if @var{AA} is upper triangular (i.e., when all the generalized eigenvalues are real, in which case the real and complex QZ coincide). @end enumerate Note: @code{qz} performs permutation balancing, but not scaling (@pxref{XREFbalance,,@code{balance}}), which may be lead to less accurate results than @code{eig}. The order of output arguments was selected for compatibility with @sc{matlab}. @seealso{eig, gsvd, balance, chol, hess, lu, qr, qzhess, schur} @end deftypefn */) { int nargin = args.length (); #if defined (DEBUG) octave_stdout << "qz: nargin = " << nargin << ", nargout = " << nargout << std::endl; #endif if (nargin < 2 || nargin > 3 || nargout > 6) print_usage (); bool complex_case = true; if (nargin == 3) { std::string opt = args(2).xstring_value ("qz: OPT must be a string"); if (opt == "real") { if (args(0).iscomplex () || args(1).iscomplex ()) { warning ("qz: ignoring 'real' option with complex matrices"); complex_case = true; } else complex_case = false; } else if (opt == "complex") complex_case = true; else error ("qz: OPT must be 'real' or 'complex'"); } #if defined (DEBUG) octave_stdout << "qz: check matrix A" << std::endl; #endif // Matrix A: check dimensions. F77_INT nn = to_f77_int (args(0).rows ()); F77_INT nc = to_f77_int (args(0).columns ()); #if defined (DEBUG) octave_stdout << "Matrix A dimensions: (" << nn << ',' << nc << ')' << std::endl; #endif if (nc != nn) err_square_matrix_required ("qz", "A"); // Matrix A: get value. Matrix aa; ComplexMatrix caa; if (args(0).iscomplex ()) caa = args(0).complex_matrix_value (); else aa = args(0).matrix_value (); #if defined (DEBUG) octave_stdout << "qz: check matrix B" << std::endl; #endif // Extract argument 2 (bb, or cbb if complex). F77_INT b_nr = to_f77_int (args(1).rows ()); F77_INT b_nc = to_f77_int (args(1).columns ()); if (nn != b_nc || nn != b_nr) ::err_nonconformant (); Matrix bb; ComplexMatrix cbb; if (args(1).iscomplex ()) cbb = args(1).complex_matrix_value (); else bb = args(1).matrix_value (); // First, declare variables used in both the real and complex cases. // FIXME: There are a lot of excess variables declared. // Probably a better way to handle this. Matrix QQ (nn, nn), ZZ (nn, nn), VR (nn, nn), VL (nn, nn); RowVector alphar (nn), alphai (nn), betar (nn); ComplexRowVector xalpha (nn), xbeta (nn); ComplexMatrix CQ (nn, nn), CZ (nn, nn), CVR (nn, nn), CVL (nn, nn); F77_INT ilo, ihi, info; char comp_q = (nargout >= 3 ? 'V' : 'N'); char comp_z = (nargout >= 4 ? 'V' : 'N'); // Initialize Q, Z to identity matrix if either is needed if (comp_q == 'V' || comp_z == 'V') { double *QQptr = QQ.rwdata (); double *ZZptr = ZZ.rwdata (); std::fill_n (QQptr, QQ.numel (), 0.0); std::fill_n (ZZptr, ZZ.numel (), 0.0); for (F77_INT i = 0; i < nn; i++) { QQ(i, i) = 1.0; ZZ(i, i) = 1.0; } } // Always perform permutation balancing. const char bal_job = 'P'; RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn); if (complex_case) { #if defined (DEBUG) if (comp_q == 'V') octave_stdout << "qz: performing balancing; CQ =\n" << CQ << std::endl; #endif if (args(0).isreal ()) caa = ComplexMatrix (aa); if (args(1).isreal ()) cbb = ComplexMatrix (bb); if (comp_q == 'V') CQ = ComplexMatrix (QQ); if (comp_z == 'V') CZ = ComplexMatrix (ZZ); F77_XFCN (zggbal, ZGGBAL, (F77_CONST_CHAR_ARG2 (&bal_job, 1), nn, F77_DBLE_CMPLX_ARG (caa.rwdata ()), nn, F77_DBLE_CMPLX_ARG (cbb.rwdata ()), nn, ilo, ihi, lscale.rwdata (), rscale.rwdata (), work.rwdata (), info F77_CHAR_ARG_LEN (1))); } else { #if defined (DEBUG) if (comp_q == 'V') octave_stdout << "qz: performing balancing; QQ =\n" << QQ << std::endl; #endif F77_XFCN (dggbal, DGGBAL, (F77_CONST_CHAR_ARG2 (&bal_job, 1), nn, aa.rwdata (), nn, bb.rwdata (), nn, ilo, ihi, lscale.rwdata (), rscale.rwdata (), work.rwdata (), info F77_CHAR_ARG_LEN (1))); } // Only permutation balance above is done. Skip scaling balance. char qz_job = (nargout < 2 ? 'E' : 'S'); if (complex_case) { // Complex case. // The QR decomposition of cbb. math::qr<ComplexMatrix> cbqr (cbb); // The R matrix of QR decomposition for cbb. cbb = cbqr.R (); // (Q*)caa for following work. caa = (cbqr.Q ().hermitian ()) * caa; CQ = CQ * cbqr.Q (); F77_XFCN (zgghrd, ZGGHRD, (F77_CONST_CHAR_ARG2 (&comp_q, 1), F77_CONST_CHAR_ARG2 (&comp_z, 1), nn, ilo, ihi, F77_DBLE_CMPLX_ARG (caa.rwdata ()), nn, F77_DBLE_CMPLX_ARG (cbb.rwdata ()), nn, F77_DBLE_CMPLX_ARG (CQ.rwdata ()), nn, F77_DBLE_CMPLX_ARG (CZ.rwdata ()), nn, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); ComplexRowVector cwork (nn); F77_XFCN (zhgeqz, ZHGEQZ, (F77_CONST_CHAR_ARG2 (&qz_job, 1), F77_CONST_CHAR_ARG2 (&comp_q, 1), F77_CONST_CHAR_ARG2 (&comp_z, 1), nn, ilo, ihi, F77_DBLE_CMPLX_ARG (caa.rwdata ()), nn, F77_DBLE_CMPLX_ARG (cbb.rwdata ()), nn, F77_DBLE_CMPLX_ARG (xalpha.rwdata ()), F77_DBLE_CMPLX_ARG (xbeta.rwdata ()), F77_DBLE_CMPLX_ARG (CQ.rwdata ()), nn, F77_DBLE_CMPLX_ARG (CZ.rwdata ()), nn, F77_DBLE_CMPLX_ARG (cwork.rwdata ()), nn, rwork.rwdata (), info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); if (comp_q == 'V') { // Left eigenvector. F77_XFCN (zggbak, ZGGBAK, (F77_CONST_CHAR_ARG2 (&bal_job, 1), F77_CONST_CHAR_ARG2 ("L", 1), nn, ilo, ihi, lscale.data (), rscale.data (), nn, F77_DBLE_CMPLX_ARG (CQ.rwdata ()), nn, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); } if (comp_z == 'V') { // Right eigenvector. F77_XFCN (zggbak, ZGGBAK, (F77_CONST_CHAR_ARG2 (&bal_job, 1), F77_CONST_CHAR_ARG2 ("R", 1), nn, ilo, ihi, lscale.data (), rscale.data (), nn, F77_DBLE_CMPLX_ARG (CZ.rwdata ()), nn, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); } } else { #if defined (DEBUG) octave_stdout << "qz: performing qr decomposition of bb" << std::endl; #endif // Compute the QR factorization of bb. math::qr<Matrix> bqr (bb); #if defined (DEBUG) octave_stdout << "qz: qr (bb) done; now performing qz decomposition" << std::endl; #endif bb = bqr.R (); #if defined (DEBUG) octave_stdout << "qz: extracted bb" << std::endl; #endif aa = (bqr.Q ()).transpose () * aa; #if defined (DEBUG) octave_stdout << "qz: updated aa " << std::endl; octave_stdout << "bqr.Q () =\n" << bqr.Q () << std::endl; if (comp_q == 'V') octave_stdout << "QQ =" << QQ << std::endl; #endif if (comp_q == 'V') QQ = QQ * bqr.Q (); #if defined (DEBUG) octave_stdout << "qz: precursors done..." << std::endl; #endif #if defined (DEBUG) octave_stdout << "qz: comp_q = " << comp_q << ", comp_z = " << comp_z << std::endl; #endif // Reduce to generalized Hessenberg form. F77_XFCN (dgghrd, DGGHRD, (F77_CONST_CHAR_ARG2 (&comp_q, 1), F77_CONST_CHAR_ARG2 (&comp_z, 1), nn, ilo, ihi, aa.rwdata (), nn, bb.rwdata (), nn, QQ.rwdata (), nn, ZZ.rwdata (), nn, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); // Check if just computing generalized eigenvalues, // or if we're actually computing the decomposition. // Reduce to generalized Schur form. F77_XFCN (dhgeqz, DHGEQZ, (F77_CONST_CHAR_ARG2 (&qz_job, 1), F77_CONST_CHAR_ARG2 (&comp_q, 1), F77_CONST_CHAR_ARG2 (&comp_z, 1), nn, ilo, ihi, aa.rwdata (), nn, bb.rwdata (), nn, alphar.rwdata (), alphai.rwdata (), betar.rwdata (), QQ.rwdata (), nn, ZZ.rwdata (), nn, work.rwdata (), nn, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); if (comp_q == 'V') { F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&bal_job, 1), F77_CONST_CHAR_ARG2 ("L", 1), nn, ilo, ihi, lscale.data (), rscale.data (), nn, QQ.rwdata (), nn, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); #if defined (DEBUG) if (comp_q == 'V') octave_stdout << "qz: balancing done; QQ=\n" << QQ << std::endl; #endif } // then right if (comp_z == 'V') { F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&bal_job, 1), F77_CONST_CHAR_ARG2 ("R", 1), nn, ilo, ihi, lscale.data (), rscale.data (), nn, ZZ.rwdata (), nn, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); #if defined (DEBUG) if (comp_z == 'V') octave_stdout << "qz: balancing done; ZZ=\n" << ZZ << std::endl; #endif } } // Right, left eigenvector matrices. if (nargout >= 5) { // Which side to compute? char side = (nargout == 5 ? 'R' : 'B'); // Compute all of them and backtransform char howmany = 'B'; // Dummy pointer; select is not used. F77_INT *select = nullptr; if (complex_case) { CVL = CQ; CVR = CZ; ComplexRowVector cwork2 (2 * nn); RowVector rwork2 (8 * nn); F77_INT m; F77_XFCN (ztgevc, ZTGEVC, (F77_CONST_CHAR_ARG2 (&side, 1), F77_CONST_CHAR_ARG2 (&howmany, 1), select, nn, F77_DBLE_CMPLX_ARG (caa.rwdata ()), nn, F77_DBLE_CMPLX_ARG (cbb.rwdata ()), nn, F77_DBLE_CMPLX_ARG (CVL.rwdata ()), nn, F77_DBLE_CMPLX_ARG (CVR.rwdata ()), nn, nn, m, F77_DBLE_CMPLX_ARG (cwork2.rwdata ()), rwork2.rwdata (), info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); } else { #if defined (DEBUG) octave_stdout << "qz: computing generalized eigenvectors" << std::endl; #endif VL = QQ; VR = ZZ; F77_INT m; F77_XFCN (dtgevc, DTGEVC, (F77_CONST_CHAR_ARG2 (&side, 1), F77_CONST_CHAR_ARG2 (&howmany, 1), select, nn, aa.rwdata (), nn, bb.rwdata (), nn, VL.rwdata (), nn, VR.rwdata (), nn, nn, m, work.rwdata (), info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); } } octave_value_list retval (nargout); switch (nargout) { case 6: // Return left eigenvectors. if (complex_case) retval(5) = CVL; else retval(5) = VL; OCTAVE_FALLTHROUGH; case 5: // Return right eigenvectors. if (complex_case) retval(4) = CVR; else retval(4) = VR; OCTAVE_FALLTHROUGH; case 4: if (complex_case) retval(3) = CZ; else retval(3) = ZZ; OCTAVE_FALLTHROUGH; case 3: if (complex_case) retval(2) = CQ.hermitian (); else retval(2) = QQ.transpose (); OCTAVE_FALLTHROUGH; case 2: case 1: case 0: { if (complex_case) { #if defined (DEBUG) octave_stdout << "qz: retval(1) = cbb =\n"; octave_print_internal (octave_stdout, cbb); octave_stdout << "\nqz: retval(0) = caa =\n"; octave_print_internal (octave_stdout, caa); octave_stdout << std::endl; #endif retval(1) = cbb; retval(0) = caa; } else { #if defined (DEBUG) octave_stdout << "qz: retval(1) = bb =\n"; octave_print_internal (octave_stdout, bb); octave_stdout << "\nqz: retval(0) = aa =\n"; octave_print_internal (octave_stdout, aa); octave_stdout << std::endl; #endif retval(1) = bb; retval(0) = aa; } } break; } // FIXME: The API for qz changed in version 9. // These warnings can be removed in Octave version 11. if (nargout == 1) { warning_with_id ("Octave:qz:single-arg-out", "qz: requesting a single output argument no longer returns eigenvalues since version 9"); disable_warning ("Octave:qz:single-arg-out"); } if (nargin == 2 && args(0).isreal () && args(1).isreal () && retval(0).iscomplex ()) { warning_with_id ("Octave:qz:complex-default", "qz: returns the complex QZ by default on real matrices since version 9"); disable_warning ("Octave:qz:complex-default"); } #if defined (DEBUG) octave_stdout << "qz: exiting (at long last)" << std::endl; #endif return retval; } /* ## Example 7.7.3 in Golub & Van Loan (3rd ed.; not present in 4th ed.) ## The generalized eigenvalues are real, so the real and complex QZ coincide. %!test %! a = [ 10 1 2; %! 1 2 -1; %! 1 1 2 ]; %! b = reshape (1:9, 3,3); %! [aa, bb, q, z, v, w] = qz (a, b); %! assert (isreal (aa) && isreal (bb) && isreal (q) && isreal (z) ... %! && isreal (v) && isreal (w)); %! assert (istriu (aa) && istriu (bb)); %! assert (q * q', eye(3), 1e-15); %! assert (z * z', eye(3), 1e-15); %! assert (q * a * z, aa, norm (aa) * 1e-14); %! assert (q * b * z, bb, norm (bb) * 1e-14); %! assert (a * v * diag (diag (bb)), b * v * diag (diag (aa)), 1e-13); %! assert (diag (diag (bb)) * w' * a, diag (diag (aa)) * w' * b, 1e-13); ## An example with real matrices where some generalized eigenvalues are ## complex, so the real and complex QZ differ. %!test %! A = [ -1.03428 0.24929 0.43205 -0.12860; %! 1.16228 0.27870 2.12954 0.69250; %! -0.51524 -0.34939 -0.77820 2.13721; %! -1.32941 2.11870 0.72005 1.00835 ]; %! B = [ 1.407302 -0.632956 -0.360628 0.068534; %! 0.149898 0.298248 0.991777 0.023652; %! 0.169281 -0.405205 -1.775834 1.511730; %! 0.717770 1.291390 -1.766607 -0.531352 ]; %! [CAA, CBB, CQ, CZ, CV, CW] = qz (A, B, 'complex'); %! assert (iscomplex (CAA) && iscomplex (CBB) && iscomplex (CQ) ... %! && iscomplex (CZ) && iscomplex (CV) && iscomplex (CW)); %! assert (istriu (CAA) && istriu (CBB)); %! assert (CQ * CQ', eye(4), 1e-14); %! assert (CZ * CZ', eye(4), 1e-14); %! assert (CQ * A * CZ, CAA, norm (CAA) * 1e-14); %! assert (CQ * B * CZ, CBB, norm (CBB) * 1e-14); %! assert (A * CV * diag (diag (CBB)), B * CV * diag (diag (CAA)), 1e-13); %! assert (diag (diag (CBB)) * CW' * A, diag (diag (CAA)) * CW' * B, 1e-13); %! [AA, BB, Q, Z, V, W] = qz (A, B, 'real'); %! assert (isreal (AA) && isreal (BB) && isreal (Q) && isreal (Z) ... %! && isreal (V) && isreal (W)); %! assert (istriu (BB)); %! assert (Q * Q', eye(4), 1e-14); %! assert (Z * Z', eye(4), 1e-14); %! assert (Q * A * Z, AA, norm (AA) * 1e-14); %! assert (Q * B * Z, BB, norm (BB) * 1e-14); ## Test input validation %!error <Invalid call> qz () %!error <Invalid call> qz (1) %!error <Invalid call> qz (1,2,3,4) %!error <Invalid call> [r1,r2,r3,r4,r5,r6,r7] = qz (1,2) %!error <OPT must be a string> qz (1,2, 3) %!error <OPT must be 'real' or 'complex'> qz (1,2, 'foobar') %!warning <ignoring 'real' option with complex matrices> qz (2i, 3, 'real'); %!warning <ignoring 'real' option with complex matrices> qz (2, 3i, 'real'); */ OCTAVE_END_NAMESPACE(octave)