qz.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 1998-2012 A. S. Hodel
00004 
00005 This file is part of Octave.
00006 
00007 Octave is free software; you can redistribute it and/or modify it
00008 under the terms of the GNU General Public License as published by the
00009 Free Software Foundation; either version 3 of the License, or (at your
00010 option) any later version.
00011 
00012 Octave is distributed in the hope that it will be useful, but WITHOUT
00013 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00014 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00015 for more details.
00016 
00017 You should have received a copy of the GNU General Public License
00018 along with Octave; see the file COPYING.  If not, see
00019 <http://www.gnu.org/licenses/>.
00020 
00021 */
00022 
00023 // Generalized eigenvalue balancing via LAPACK
00024 
00025 // Author: A. S. Hodel <scotte@eng.auburn.edu>
00026 
00027 #undef DEBUG
00028 #undef DEBUG_SORT
00029 #undef DEBUG_EIG
00030 
00031 #ifdef HAVE_CONFIG_H
00032 #include <config.h>
00033 #endif
00034 
00035 #include <cfloat>
00036 
00037 #include <iostream>
00038 #include <iomanip>
00039 
00040 #include "CmplxQRP.h"
00041 #include "CmplxQR.h"
00042 #include "dbleQR.h"
00043 #include "f77-fcn.h"
00044 #include "lo-math.h"
00045 #include "quit.h"
00046 
00047 #include "defun-dld.h"
00048 #include "error.h"
00049 #include "gripes.h"
00050 #include "oct-obj.h"
00051 #include "oct-map.h"
00052 #include "ov.h"
00053 #include "pager.h"
00054 #if defined (DEBUG) || defined (DEBUG_SORT)
00055 #include "pr-output.h"
00056 #endif
00057 #include "symtab.h"
00058 #include "utils.h"
00059 #include "variables.h"
00060 
00061 typedef octave_idx_type (*sort_function) (const octave_idx_type& LSIZE, const double& ALPHA,
00062                               const double& BETA, const double& S,
00063                               const double& P);
00064 
00065 extern "C"
00066 {
00067   F77_RET_T
00068   F77_FUNC (dggbal, DGGBAL) (F77_CONST_CHAR_ARG_DECL,
00069                              const octave_idx_type& N, double* A,
00070                              const octave_idx_type& LDA, double* B,
00071                              const octave_idx_type& LDB, octave_idx_type& ILO,
00072                              octave_idx_type& IHI, double* LSCALE,
00073                              double* RSCALE, double* WORK,
00074                              octave_idx_type& INFO
00075                              F77_CHAR_ARG_LEN_DECL);
00076 
00077 F77_RET_T
00078   F77_FUNC (zggbal, ZGGBAL) (F77_CONST_CHAR_ARG_DECL,
00079                              const octave_idx_type& N, Complex* A,
00080                              const octave_idx_type& LDA, Complex* B,
00081                              const octave_idx_type& LDB, octave_idx_type& ILO,
00082                              octave_idx_type& IHI, double* LSCALE,
00083                              double* RSCALE, double* WORK,
00084                              octave_idx_type& INFO
00085                              F77_CHAR_ARG_LEN_DECL);
00086 
00087   F77_RET_T
00088   F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL,
00089                              F77_CONST_CHAR_ARG_DECL,
00090                              const octave_idx_type& N,
00091                              const octave_idx_type& ILO,
00092                              const octave_idx_type& IHI,
00093                              const double* LSCALE, const double* RSCALE,
00094                              octave_idx_type& M, double* V,
00095                              const octave_idx_type& LDV, octave_idx_type& INFO
00096                              F77_CHAR_ARG_LEN_DECL
00097                              F77_CHAR_ARG_LEN_DECL);
00098 
00099 F77_RET_T
00100   F77_FUNC (zggbak, ZGGBAK) (F77_CONST_CHAR_ARG_DECL,
00101                              F77_CONST_CHAR_ARG_DECL,
00102                              const octave_idx_type& N,
00103                              const octave_idx_type& ILO,
00104                              const octave_idx_type& IHI,
00105                              const double* LSCALE, const double* RSCALE,
00106                              octave_idx_type& M, Complex* V,
00107                              const octave_idx_type& LDV, octave_idx_type& INFO
00108                              F77_CHAR_ARG_LEN_DECL
00109                              F77_CHAR_ARG_LEN_DECL);
00110 
00111   F77_RET_T
00112   F77_FUNC (dgghrd, DGGHRD) (F77_CONST_CHAR_ARG_DECL,
00113                              F77_CONST_CHAR_ARG_DECL,
00114                              const octave_idx_type& N,
00115                              const octave_idx_type& ILO,
00116                              const octave_idx_type& IHI, double* A,
00117                              const octave_idx_type& LDA, double* B,
00118                              const octave_idx_type& LDB, double* Q,
00119                              const octave_idx_type& LDQ, double* Z,
00120                              const octave_idx_type& LDZ, octave_idx_type& INFO
00121                              F77_CHAR_ARG_LEN_DECL
00122                              F77_CHAR_ARG_LEN_DECL);
00123 
00124  F77_RET_T
00125   F77_FUNC (zgghrd, ZGGHRD) (F77_CONST_CHAR_ARG_DECL,
00126                              F77_CONST_CHAR_ARG_DECL,
00127                              const octave_idx_type& N,
00128                              const octave_idx_type& ILO,
00129                              const octave_idx_type& IHI, Complex* A,
00130                              const octave_idx_type& LDA, Complex* B,
00131                              const octave_idx_type& LDB, Complex* Q,
00132                              const octave_idx_type& LDQ, Complex* Z,
00133                              const octave_idx_type& LDZ, octave_idx_type& INFO
00134                              F77_CHAR_ARG_LEN_DECL
00135                              F77_CHAR_ARG_LEN_DECL);
00136 
00137   F77_RET_T
00138   F77_FUNC (dhgeqz, DHGEQZ) (F77_CONST_CHAR_ARG_DECL,
00139                              F77_CONST_CHAR_ARG_DECL,
00140                              F77_CONST_CHAR_ARG_DECL,
00141                              const octave_idx_type& N,
00142                              const octave_idx_type& ILO,
00143                              const octave_idx_type& IHI,
00144                              double* A, const octave_idx_type& LDA, double* B,
00145                              const octave_idx_type& LDB, double* ALPHAR,
00146                              double* ALPHAI, double* BETA, double* Q,
00147                              const octave_idx_type& LDQ, double* Z,
00148                              const octave_idx_type& LDZ, double* WORK,
00149                              const octave_idx_type& LWORK,
00150                              octave_idx_type& INFO
00151                              F77_CHAR_ARG_LEN_DECL
00152                              F77_CHAR_ARG_LEN_DECL
00153                              F77_CHAR_ARG_LEN_DECL);
00154 
00155 F77_RET_T
00156   F77_FUNC (zhgeqz, ZHGEQZ) (F77_CONST_CHAR_ARG_DECL,
00157                              F77_CONST_CHAR_ARG_DECL,
00158                              F77_CONST_CHAR_ARG_DECL,
00159                              const octave_idx_type& N,
00160                              const octave_idx_type& ILO,
00161                              const octave_idx_type& IHI,
00162                              Complex* A, const octave_idx_type& LDA,
00163                              Complex* B, const octave_idx_type& LDB,
00164                              Complex* ALPHA, Complex* BETA, Complex* CQ,
00165                              const octave_idx_type& LDQ,
00166                              Complex* CZ, const octave_idx_type& LDZ,
00167                              Complex* WORK, const octave_idx_type& LWORK,
00168                              double* RWORK, octave_idx_type& INFO
00169                              F77_CHAR_ARG_LEN_DECL
00170                              F77_CHAR_ARG_LEN_DECL
00171                              F77_CHAR_ARG_LEN_DECL);
00172 
00173   F77_RET_T
00174   F77_FUNC (dlag2, DLAG2) (const double* A, const octave_idx_type& LDA,
00175                            const double* B, const octave_idx_type& LDB,
00176                            const double& SAFMIN, double& SCALE1,
00177                            double& SCALE2, double& WR1, double& WR2,
00178                            double& WI);
00179 
00180   // Van Dooren's code (netlib.org: toms/590) for reordering
00181   // GEP.  Only processes Z, not Q.
00182   F77_RET_T
00183   F77_FUNC (dsubsp, DSUBSP) (const octave_idx_type& NMAX,
00184                              const octave_idx_type& N, double* A,
00185                              double* B, double* Z, sort_function,
00186                              const double& EPS, octave_idx_type& NDIM,
00187                              octave_idx_type& FAIL, octave_idx_type* IND);
00188 
00189   // Documentation for DTGEVC incorrectly states that VR, VL are
00190   // complex*16; they are declared in DTGEVC as double precision
00191   // (probably a cut and paste problem fro ZTGEVC).
00192   F77_RET_T
00193   F77_FUNC (dtgevc, DTGEVC) (F77_CONST_CHAR_ARG_DECL,
00194                              F77_CONST_CHAR_ARG_DECL,
00195                              octave_idx_type* SELECT,
00196                              const octave_idx_type& N, double* A,
00197                              const octave_idx_type& LDA, double* B,
00198                              const octave_idx_type& LDB, double* VL,
00199                              const octave_idx_type& LDVL, double* VR,
00200                              const octave_idx_type& LDVR,
00201                              const octave_idx_type& MM, octave_idx_type& M,
00202                              double* WORK, octave_idx_type& INFO
00203                              F77_CHAR_ARG_LEN_DECL
00204                              F77_CHAR_ARG_LEN_DECL);
00205 
00206 F77_RET_T
00207   F77_FUNC (ztgevc, ZTGEVC) (F77_CONST_CHAR_ARG_DECL,
00208                              F77_CONST_CHAR_ARG_DECL,
00209                              octave_idx_type* SELECT,
00210                              const octave_idx_type& N, const Complex* A,
00211                              const octave_idx_type& LDA,const Complex* B,
00212                              const octave_idx_type& LDB, Complex* xVL,
00213                              const octave_idx_type& LDVL, Complex* xVR,
00214                              const octave_idx_type& LDVR,
00215                              const octave_idx_type& MM, octave_idx_type& M,
00216                              Complex* CWORK, double* RWORK,
00217                              octave_idx_type& INFO
00218                              F77_CHAR_ARG_LEN_DECL
00219                              F77_CHAR_ARG_LEN_DECL);
00220 
00221   F77_RET_T
00222   F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG_DECL,
00223                                double& retval
00224                                F77_CHAR_ARG_LEN_DECL);
00225 
00226   F77_RET_T
00227   F77_FUNC (xdlange, XDLANGE) (F77_CONST_CHAR_ARG_DECL,
00228                                const octave_idx_type&,
00229                                const octave_idx_type&, const double*,
00230                                const octave_idx_type&, double*, double&
00231                                F77_CHAR_ARG_LEN_DECL);
00232 }
00233 
00234 // fcrhp, fin, fout, folhp:
00235 // routines for ordering of generalized eigenvalues
00236 // return 1 if  test is passed, 0 otherwise
00237 //    fin: |lambda| < 1
00238 //    fout: |lambda| >= 1
00239 //    fcrhp: real(lambda) >= 0
00240 //    folhp: real(lambda) < 0
00241 
00242 static octave_idx_type
00243 fcrhp (const octave_idx_type& lsize, const double& alpha,
00244        const double& beta, const double& s, const double&)
00245 {
00246   if (lsize == 1)
00247     return (alpha * beta >= 0 ? 1 : -1);
00248   else
00249     return (s >= 0 ? 1 : -1);
00250 }
00251 
00252 static octave_idx_type
00253 fin (const octave_idx_type& lsize, const double& alpha,
00254      const double& beta, const double&, const double& p)
00255 {
00256   octave_idx_type retval;
00257 
00258   if (lsize == 1)
00259     retval = (fabs (alpha) < fabs (beta) ? 1 : -1);
00260   else
00261     retval = (fabs (p) < 1 ? 1 : -1);
00262 
00263 #ifdef DEBUG
00264   std::cout << "qz: fin: retval=" << retval << std::endl;
00265 #endif
00266 
00267   return retval;
00268 }
00269 
00270 static octave_idx_type
00271 folhp (const octave_idx_type& lsize, const double& alpha,
00272        const double& beta, const double& s, const double&)
00273 {
00274   if (lsize == 1)
00275     return (alpha * beta < 0 ? 1 : -1);
00276   else
00277     return (s < 0 ? 1 : -1);
00278 }
00279 
00280 static octave_idx_type
00281 fout (const octave_idx_type& lsize, const double& alpha,
00282       const double& beta, const double&, const double& p)
00283 {
00284   if (lsize == 1)
00285     return (fabs (alpha) >= fabs (beta) ? 1 : -1);
00286   else
00287     return (fabs (p) >= 1 ? 1 : -1);
00288 }
00289 
00290 
00291 //FIXME: Matlab does not produce lambda as the first output argument.
00292 //       Compatibility problem?
00293 DEFUN_DLD (qz, args, nargout,
00294   "-*- texinfo -*-\n\
00295 @deftypefn  {Loadable Function} {@var{lambda} =} qz (@var{A}, @var{B})\n\
00296 @deftypefnx {Loadable Function} {@var{lambda} =} qz (@var{A}, @var{B}, @var{opt})\n\
00297 QZ@tie{}decomposition of the generalized eigenvalue problem\n\
00298 (@math{A x = s B x}).  There are three ways to call this function:\n\
00299 @enumerate\n\
00300 @item @code{@var{lambda} = qz (@var{A}, @var{B})}\n\
00301 \n\
00302 Computes the generalized eigenvalues\n\
00303 @tex\n\
00304 $\\lambda$\n\
00305 @end tex\n\
00306 @ifnottex\n\
00307 @var{lambda}\n\
00308 @end ifnottex\n\
00309 of @math{(A - s B)}.\n\
00310 \n\
00311 @item @code{[AA, BB, Q, Z, V, W, @var{lambda}] = qz (@var{A}, @var{B})}\n\
00312 \n\
00313 Computes QZ@tie{}decomposition, generalized eigenvectors, and\n\
00314 generalized eigenvalues of @math{(A - s B)}\n\
00315 @tex\n\
00316 $$ AV = BV{ \\rm diag }(\\lambda) $$\n\
00317 $$ W^T A = { \\rm diag }(\\lambda)W^T B $$\n\
00318 $$ AA = Q^T AZ, BB = Q^T BZ $$\n\
00319 @end tex\n\
00320 @ifnottex\n\
00321 \n\
00322 @example\n\
00323 @group\n\
00324 \n\
00325     A * V = B * V * diag (@var{lambda})\n\
00326     W' * A = diag (@var{lambda}) * W' * B\n\
00327     AA = Q * A * Z, BB = Q * B * Z\n\
00328 \n\
00329 @end group\n\
00330 @end example\n\
00331 \n\
00332 @end ifnottex\n\
00333 with @var{Q} and @var{Z} orthogonal (unitary)= @var{I}\n\
00334 \n\
00335 @item @code{[AA,BB,Z@{, @var{lambda}@}] = qz (@var{A}, @var{B}, @var{opt})}\n\
00336 \n\
00337 As in form [2], but allows ordering of generalized eigenpairs\n\
00338 for (e.g.) solution of discrete time algebraic Riccati equations.\n\
00339 Form 3 is not available for complex matrices, and does not compute\n\
00340 the generalized eigenvectors @var{V}, @var{W}, nor the orthogonal matrix\n\
00341 @var{Q}.\n\
00342 \n\
00343 @table @var\n\
00344 @item opt\n\
00345 for ordering eigenvalues of the GEP pencil.  The leading block\n\
00346 of the revised pencil contains all eigenvalues that satisfy:\n\
00347 @table @asis\n\
00348 @item \"N\"\n\
00349 = unordered (default)\n\
00350 \n\
00351 @item \"S\"\n\
00352 = small: leading block has all |lambda| @leq{} 1\n\
00353 \n\
00354 @item \"B\"\n\
00355 = big: leading block has all |lambda| @geq{} 1\n\
00356 \n\
00357 @item \"-\"\n\
00358 = negative real part: leading block has all eigenvalues\n\
00359 in the open left half-plane\n\
00360 \n\
00361 @item \"+\"\n\
00362 = non-negative real part: leading block has all eigenvalues\n\
00363 in the closed right half-plane\n\
00364 @end table\n\
00365 @end table\n\
00366 @end enumerate\n\
00367 \n\
00368 Note: @code{qz} performs permutation balancing, but not scaling\n\
00369 (@pxref{doc-balance}).  The order of output arguments was selected for\n\
00370 compatibility with @sc{matlab}.\n\
00371 @seealso{balance, eig, schur}\n\
00372 @end deftypefn")
00373 {
00374   octave_value_list retval;
00375   int nargin = args.length ();
00376 
00377 #ifdef DEBUG
00378   std::cout << "qz: nargin = " << nargin << ", nargout = " << nargout << std::endl;
00379 #endif
00380 
00381   if (nargin < 2 || nargin > 3 || nargout > 7)
00382     {
00383       print_usage ();
00384       return retval;
00385     }
00386   else if (nargin == 3 && (nargout < 3 || nargout > 4))
00387     {
00388       error ("qz: invalid number of output arguments for form [3] call");
00389       return retval;
00390     }
00391 
00392 #ifdef DEBUG
00393   std::cout << "qz: determine ordering option" << std::endl;
00394 #endif
00395 
00396   // Determine ordering option.
00397   volatile char ord_job = 0;
00398   static double safmin;
00399 
00400   if (nargin == 2)
00401     ord_job = 'N';
00402   else if (!args(2).is_string ())
00403     {
00404       error ("qz: OPT must be a string");
00405       return retval;
00406     }
00407   else
00408     {
00409       std::string tmp = args(2).string_value ();
00410 
00411       if (! tmp.empty ())
00412         ord_job = tmp[0];
00413 
00414       if (! (ord_job == 'N' || ord_job == 'n'
00415              || ord_job == 'S' || ord_job == 's'
00416              || ord_job == 'B' || ord_job == 'b'
00417              || ord_job == '+' || ord_job == '-'))
00418         {
00419           error ("qz: invalid order option");
00420           return retval;
00421         }
00422 
00423       // overflow constant required by dlag2
00424       F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("S", 1),
00425                                    safmin
00426                                    F77_CHAR_ARG_LEN (1));
00427 
00428 #ifdef DEBUG_EIG
00429       std::cout << "qz: initial value of safmin=" << setiosflags (std::ios::scientific)
00430            << safmin << std::endl;
00431 #endif
00432 
00433       // Some machines (e.g., DEC alpha) get safmin = 0;
00434       // for these, use eps instead to avoid problems in dlag2.
00435       if (safmin == 0)
00436         {
00437 #ifdef DEBUG_EIG
00438           std::cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << std::endl;
00439 #endif
00440 
00441           F77_FUNC (xdlamch, XDLAMCH) (F77_CONST_CHAR_ARG2 ("E", 1),
00442                                        safmin
00443                                        F77_CHAR_ARG_LEN (1));
00444 
00445 #ifdef DEBUG_EIG
00446           std::cout << "qz: safmin set to " << setiosflags (std::ios::scientific)
00447                << safmin << std::endl;
00448 #endif
00449         }
00450     }
00451 
00452 #ifdef DEBUG
00453   std::cout << "qz: check argument 1" << std::endl;
00454 #endif
00455 
00456   // Argument 1: check if it's o.k. dimensioned.
00457   octave_idx_type nn = args(0).rows ();
00458 
00459 #ifdef DEBUG
00460   std::cout << "argument 1 dimensions: (" << nn << "," << args(0).columns () << ")"
00461        << std::endl;
00462 #endif
00463 
00464   int arg_is_empty = empty_arg ("qz", nn, args(0).columns ());
00465 
00466   if (arg_is_empty < 0)
00467     {
00468       gripe_empty_arg ("qz: parameter 1", 0);
00469       return retval;
00470     }
00471   else if (arg_is_empty > 0)
00472     {
00473       gripe_empty_arg ("qz: parameter 1; continuing", 0);
00474       return octave_value_list (2, Matrix ());
00475     }
00476   else if (args(0).columns () != nn)
00477     {
00478       gripe_square_matrix_required ("qz");
00479       return retval;
00480     }
00481 
00482   // Argument 1: dimensions look good; get the value.
00483   Matrix aa;
00484   ComplexMatrix caa;
00485 
00486   if (args(0).is_complex_type ())
00487     caa = args(0).complex_matrix_value ();
00488   else
00489     aa = args(0).matrix_value ();
00490 
00491   if (error_state)
00492     return retval;
00493 
00494 #ifdef DEBUG
00495   std::cout << "qz: check argument 2" << std::endl;
00496 #endif
00497 
00498   // Extract argument 2 (bb, or cbb if complex).
00499   if ((nn != args(1).columns ()) || (nn != args(1).rows ()))
00500     {
00501       gripe_nonconformant ();
00502       return retval;
00503     }
00504 
00505   Matrix bb;
00506   ComplexMatrix cbb;
00507 
00508   if (args(1).is_complex_type ())
00509     cbb = args(1).complex_matrix_value ();
00510   else
00511     bb = args(1).matrix_value ();
00512 
00513   if (error_state)
00514     return retval;
00515 
00516   // Both matrices loaded, now let's check what kind of arithmetic:
00517   // declared volatile to avoid compiler warnings about long jumps,
00518   // vforks.
00519 
00520   volatile int complex_case
00521     = (args(0).is_complex_type () || args(1).is_complex_type ());
00522 
00523   if (nargin == 3 && complex_case)
00524     {
00525       error ("qz: cannot re-order complex qz decomposition");
00526       return retval;
00527     }
00528 
00529   // First, declare variables used in both the real and complex case.
00530   Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn);
00531   RowVector alphar(nn), alphai(nn), betar(nn);
00532   ComplexRowVector xalpha(nn), xbeta(nn);
00533   ComplexMatrix CQ(nn,nn), CZ(nn,nn), CVR(nn,nn), CVL(nn,nn);
00534   octave_idx_type ilo, ihi, info;
00535   char compq = (nargout >= 3 ? 'V' : 'N');
00536   char compz = (nargout >= 4 ? 'V' : 'N');
00537 
00538   // Initialize Q, Z to identity if we need either of them.
00539   if (compq == 'V' || compz == 'V')
00540     for (octave_idx_type ii = 0; ii < nn; ii++)
00541       for (octave_idx_type jj = 0; jj < nn; jj++)
00542         {
00543           OCTAVE_QUIT;
00544           QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0);
00545         }
00546 
00547   // Always perform permutation balancing.
00548   const char bal_job = 'P';
00549   RowVector lscale (nn), rscale (nn), work (6 * nn), rwork (nn);
00550 
00551   if (complex_case)
00552     {
00553 #ifdef DEBUG
00554       if (compq == 'V')
00555         std::cout << "qz: performing balancing; CQ=" << std::endl << CQ << std::endl;
00556 #endif
00557       if (args(0).is_real_type ())
00558         caa = ComplexMatrix (aa);
00559 
00560       if (args(1).is_real_type ())
00561         cbb = ComplexMatrix (bb);
00562 
00563       if (compq == 'V')
00564         CQ = ComplexMatrix (QQ);
00565 
00566       if (compz == 'V')
00567         CZ = ComplexMatrix (ZZ);
00568 
00569       F77_XFCN (zggbal, ZGGBAL,
00570                 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
00571                  nn, caa.fortran_vec (), nn, cbb.fortran_vec (),
00572                  nn, ilo, ihi, lscale.fortran_vec (),
00573                  rscale.fortran_vec (), work.fortran_vec (), info
00574                  F77_CHAR_ARG_LEN (1)));
00575     }
00576   else
00577     {
00578 #ifdef DEBUG
00579       if (compq == 'V')
00580         std::cout << "qz: performing balancing; QQ=" << std::endl << QQ << std::endl;
00581 #endif
00582 
00583       F77_XFCN (dggbal, DGGBAL,
00584                 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
00585                  nn, aa.fortran_vec (), nn, bb.fortran_vec (),
00586                  nn, ilo, ihi, lscale.fortran_vec (),
00587                  rscale.fortran_vec (), work.fortran_vec (), info
00588                  F77_CHAR_ARG_LEN (1)));
00589     }
00590 
00591   // Since we just want the balancing matrices, we can use dggbal
00592   // for both the real and complex cases; left first
00593 
00594 #if 0
00595   if (compq == 'V')
00596     {
00597       F77_XFCN (dggbak, DGGBAK,
00598                 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
00599                  F77_CONST_CHAR_ARG2 ("L", 1),
00600                  nn, ilo, ihi, lscale.data (), rscale.data (),
00601                  nn, QQ.fortran_vec (), nn, info
00602                  F77_CHAR_ARG_LEN (1)
00603                  F77_CHAR_ARG_LEN (1)));
00604 
00605 #ifdef DEBUG
00606       if (compq == 'V')
00607         std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl;
00608 #endif
00609   }
00610 
00611   // then right
00612   if (compz == 'V')
00613     {
00614       F77_XFCN (dggbak, DGGBAK,
00615                 (F77_CONST_CHAR_ARG2 (&bal_job, 1),
00616                  F77_CONST_CHAR_ARG2 ("R", 1),
00617                  nn, ilo, ihi, lscale.data (), rscale.data (),
00618                  nn, ZZ.fortran_vec (), nn, info
00619                  F77_CHAR_ARG_LEN (1)
00620                  F77_CHAR_ARG_LEN (1)));
00621 
00622 #ifdef DEBUG
00623       if (compz == 'V')
00624         std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl;
00625 #endif
00626     }
00627 #endif
00628 
00629   static char qz_job;
00630   qz_job = (nargout < 2 ? 'E' : 'S');
00631 
00632   if (complex_case)
00633     {
00634       // Complex case.
00635 
00636       // The QR decomposition of cbb.
00637       ComplexQR cbqr (cbb);
00638       // The R matrix of QR decomposition for cbb.
00639       cbb = cbqr.R ();
00640       // (Q*)caa for following work.
00641       caa = (cbqr.Q ().hermitian ()) * caa;
00642       CQ = CQ * cbqr.Q ();
00643 
00644       F77_XFCN (zgghrd, ZGGHRD,
00645                 (F77_CONST_CHAR_ARG2 (&compq, 1),
00646                  F77_CONST_CHAR_ARG2 (&compz, 1),
00647                  nn, ilo, ihi, caa.fortran_vec (),
00648                  nn, cbb.fortran_vec (), nn, CQ.fortran_vec (), nn,
00649                  CZ.fortran_vec (), nn, info
00650                  F77_CHAR_ARG_LEN (1)
00651                  F77_CHAR_ARG_LEN (1)));
00652 
00653       ComplexRowVector cwork (1 * nn);
00654 
00655       F77_XFCN (zhgeqz, ZHGEQZ,
00656                 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
00657                  F77_CONST_CHAR_ARG2 (&compq, 1),
00658                  F77_CONST_CHAR_ARG2 (&compz, 1),
00659                  nn, ilo, ihi,
00660                  caa.fortran_vec (), nn,
00661                  cbb.fortran_vec (),nn,
00662                  xalpha.fortran_vec (), xbeta.fortran_vec (),
00663                  CQ.fortran_vec (), nn,
00664                  CZ.fortran_vec (), nn,
00665                  cwork.fortran_vec (), nn, rwork.fortran_vec (), info
00666                  F77_CHAR_ARG_LEN (1)
00667                  F77_CHAR_ARG_LEN (1)
00668                  F77_CHAR_ARG_LEN (1)));
00669 
00670       if (compq == 'V')
00671         {
00672           // Left eigenvector.
00673           F77_XFCN (zggbak, ZGGBAK,
00674                     (F77_CONST_CHAR_ARG2 (&bal_job, 1),
00675                      F77_CONST_CHAR_ARG2 ("L", 1),
00676                      nn, ilo, ihi, lscale.data (), rscale.data (),
00677                      nn, CQ.fortran_vec (), nn, info
00678                      F77_CHAR_ARG_LEN (1)
00679                      F77_CHAR_ARG_LEN (1)));
00680         }
00681 
00682       // Right eigenvector.
00683       if (compz == 'V')
00684         {
00685           F77_XFCN (zggbak, ZGGBAK,
00686                     (F77_CONST_CHAR_ARG2 (&bal_job, 1),
00687                      F77_CONST_CHAR_ARG2 ("R", 1),
00688                      nn, ilo, ihi, lscale.data (), rscale.data (),
00689                      nn, CZ.fortran_vec (), nn, info
00690                      F77_CHAR_ARG_LEN (1)
00691                      F77_CHAR_ARG_LEN (1)));
00692         }
00693 
00694     }
00695   else
00696     {
00697 #ifdef DEBUG
00698       std::cout << "qz: peforming qr decomposition of bb" << std::endl;
00699 #endif
00700 
00701       // Compute the QR factorization of bb.
00702       QR bqr (bb);
00703 
00704 #ifdef DEBUG
00705       std::cout << "qz: qr (bb) done; now peforming qz decomposition" << std::endl;
00706 #endif
00707 
00708       bb = bqr.R ();
00709 
00710 #ifdef DEBUG
00711       std::cout << "qz: extracted bb" << std::endl;
00712 #endif
00713 
00714       aa = (bqr.Q ()).transpose () * aa;
00715 
00716 #ifdef DEBUG
00717       std::cout << "qz: updated aa " << std::endl;
00718       std::cout << "bqr.Q () = " << std::endl << bqr.Q () << std::endl;
00719 
00720       if (compq == 'V')
00721         std::cout << "QQ =" << QQ << std::endl;
00722 #endif
00723 
00724       if (compq == 'V')
00725         QQ = QQ * bqr.Q ();
00726 
00727 #ifdef DEBUG
00728       std::cout << "qz: precursors done..." << std::endl;
00729 #endif
00730 
00731 #ifdef DEBUG
00732       std::cout << "qz: compq = " << compq << ", compz = " << compz << std::endl;
00733 #endif
00734 
00735       // Reduce  to generalized hessenberg form.
00736       F77_XFCN (dgghrd, DGGHRD,
00737                 (F77_CONST_CHAR_ARG2 (&compq, 1),
00738                  F77_CONST_CHAR_ARG2 (&compz, 1),
00739                  nn, ilo, ihi, aa.fortran_vec (),
00740                  nn, bb.fortran_vec (), nn, QQ.fortran_vec (), nn,
00741                  ZZ.fortran_vec (), nn, info
00742                  F77_CHAR_ARG_LEN (1)
00743                  F77_CHAR_ARG_LEN (1)));
00744 
00745       // Check if just computing generalized eigenvalues or if we're
00746       // actually computing the decomposition.
00747 
00748       // Reduce to generalized Schur form.
00749       F77_XFCN (dhgeqz, DHGEQZ,
00750                 (F77_CONST_CHAR_ARG2 (&qz_job, 1),
00751                  F77_CONST_CHAR_ARG2 (&compq, 1),
00752                  F77_CONST_CHAR_ARG2 (&compz, 1),
00753                  nn, ilo, ihi, aa.fortran_vec (), nn, bb.fortran_vec (),
00754                  nn, alphar.fortran_vec (), alphai.fortran_vec (),
00755                  betar.fortran_vec (), QQ.fortran_vec (), nn,
00756                  ZZ.fortran_vec (), nn, work.fortran_vec (), nn, info
00757                  F77_CHAR_ARG_LEN (1)
00758                  F77_CHAR_ARG_LEN (1)
00759                  F77_CHAR_ARG_LEN (1)));
00760 
00761       if (compq == 'V')
00762         {
00763           F77_XFCN (dggbak, DGGBAK,
00764                     (F77_CONST_CHAR_ARG2 (&bal_job, 1),
00765                      F77_CONST_CHAR_ARG2 ("L", 1),
00766                      nn, ilo, ihi, lscale.data (), rscale.data (),
00767                      nn, QQ.fortran_vec (), nn, info
00768                      F77_CHAR_ARG_LEN (1)
00769                      F77_CHAR_ARG_LEN (1)));
00770 
00771 #ifdef DEBUG
00772           if (compq == 'V')
00773             std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl;
00774 #endif
00775         }
00776 
00777   // then right
00778       if (compz == 'V')
00779         {
00780            F77_XFCN (dggbak, DGGBAK,
00781                      (F77_CONST_CHAR_ARG2 (&bal_job, 1),
00782                       F77_CONST_CHAR_ARG2 ("R", 1),
00783                       nn, ilo, ihi, lscale.data (), rscale.data (),
00784                       nn, ZZ.fortran_vec (), nn, info
00785                       F77_CHAR_ARG_LEN (1)
00786                       F77_CHAR_ARG_LEN (1)));
00787 
00788 #ifdef DEBUG
00789            if (compz == 'V')
00790              std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl;
00791 #endif
00792         }
00793 
00794     }
00795 
00796   // Order the QZ decomposition?
00797   if (! (ord_job == 'N' || ord_job == 'n'))
00798     {
00799       if (complex_case)
00800         {
00801           // Probably not needed, but better be safe.
00802           error ("qz: cannot re-order complex qz decomposition");
00803           return retval;
00804         }
00805       else
00806         {
00807 #ifdef DEBUG_SORT
00808           std::cout << "qz: ordering eigenvalues: ord_job = "
00809                     << ord_job << std::endl;
00810 #endif
00811 
00812           // Declared static to avoid vfork/long jump compiler complaints.
00813           static sort_function sort_test;
00814           sort_test = 0;
00815 
00816           switch (ord_job)
00817             {
00818             case 'S':
00819             case 's':
00820               sort_test = &fin;
00821               break;
00822 
00823             case 'B':
00824             case 'b':
00825               sort_test = &fout;
00826               break;
00827 
00828             case '+':
00829               sort_test = &fcrhp;
00830               break;
00831 
00832             case '-':
00833               sort_test = &folhp;
00834               break;
00835 
00836             default:
00837               // Invalid order option (should never happen, since we
00838               // checked the options at the top).
00839               panic_impossible ();
00840               break;
00841             }
00842 
00843           octave_idx_type ndim, fail;
00844           double inf_norm;
00845 
00846           F77_XFCN (xdlange, XDLANGE,
00847                     (F77_CONST_CHAR_ARG2 ("I", 1),
00848                      nn, nn, aa.data (), nn, work.fortran_vec (), inf_norm
00849                      F77_CHAR_ARG_LEN (1)));
00850 
00851           double eps = DBL_EPSILON * inf_norm * nn;
00852 
00853 #ifdef DEBUG_SORT
00854           std::cout << "qz: calling dsubsp: aa=" << std::endl;
00855           octave_print_internal (std::cout, aa, 0);
00856           std::cout << std::endl << "bb="  << std::endl;
00857           octave_print_internal (std::cout, bb, 0);
00858           if (compz == 'V')
00859             {
00860               std::cout << std::endl << "ZZ="  << std::endl;
00861               octave_print_internal (std::cout, ZZ, 0);
00862             }
00863           std::cout << std::endl;
00864           std::cout << "alphar = " << std::endl;
00865           octave_print_internal (std::cout, (Matrix) alphar, 0);
00866           std::cout << std::endl << "alphai = " << std::endl;
00867           octave_print_internal (std::cout, (Matrix) alphai, 0);
00868           std::cout << std::endl << "beta = " << std::endl;
00869           octave_print_internal (std::cout, (Matrix) betar, 0);
00870           std::cout << std::endl;
00871 #endif
00872 
00873           Array<octave_idx_type> ind (dim_vector (nn, 1));
00874 
00875           F77_XFCN (dsubsp, DSUBSP,
00876                     (nn, nn, aa.fortran_vec (), bb.fortran_vec (),
00877                      ZZ.fortran_vec (), sort_test, eps, ndim, fail,
00878                      ind.fortran_vec ()));
00879 
00880 #ifdef DEBUG
00881           std::cout << "qz: back from dsubsp: aa=" << std::endl;
00882           octave_print_internal (std::cout, aa, 0);
00883           std::cout << std::endl << "bb="  << std::endl;
00884           octave_print_internal (std::cout, bb, 0);
00885           if (compz == 'V')
00886             {
00887               std::cout << std::endl << "ZZ="  << std::endl;
00888               octave_print_internal (std::cout, ZZ, 0);
00889             }
00890           std::cout << std::endl;
00891 #endif
00892 
00893           // Manually update alphar, alphai, betar.
00894           static int jj;
00895 
00896           jj = 0;
00897           while (jj < nn)
00898             {
00899 #ifdef DEBUG_EIG
00900               std::cout << "computing gen eig #" << jj << std::endl;
00901 #endif
00902 
00903               // Number of zeros in this block.
00904               static int zcnt;
00905 
00906               if (jj == (nn-1))
00907                 zcnt = 1;
00908               else if (aa(jj+1,jj) == 0)
00909                 zcnt = 1;
00910               else zcnt = 2;
00911 
00912               if (zcnt == 1)
00913                 {
00914                   // Real zero.
00915 #ifdef DEBUG_EIG
00916                   std::cout << "  single gen eig:" << std::endl;
00917                   std::cout << "  alphar(" << jj << ") = " << aa(jj,jj) << std::endl;
00918                   std::cout << "  betar( " << jj << ") = " << bb(jj,jj) << std::endl;
00919                   std::cout << "  alphai(" << jj << ") = 0" << std::endl;
00920 #endif
00921 
00922                   alphar(jj) = aa(jj,jj);
00923                   alphai(jj) = 0;
00924                   betar(jj) = bb(jj,jj);
00925                 }
00926               else
00927                 {
00928                   // Complex conjugate pair.
00929 #ifdef DEBUG_EIG
00930                   std::cout << "qz: calling dlag2:" << std::endl;
00931                   std::cout << "safmin="
00932                        << setiosflags (std::ios::scientific) << safmin << std::endl;
00933 
00934                   for (int idr = jj; idr <= jj+1; idr++)
00935                     {
00936                       for (int idc = jj; idc <= jj+1; idc++)
00937                         {
00938                           std::cout << "aa(" << idr << "," << idc << ")="
00939                                << aa(idr,idc) << std::endl;
00940                           std::cout << "bb(" << idr << "," << idc << ")="
00941                                << bb(idr,idc) << std::endl;
00942                         }
00943                     }
00944 #endif
00945 
00946                   // FIXME -- probably should be using
00947                   // fortran_vec instead of &aa(jj,jj) here.
00948 
00949                   double scale1, scale2, wr1, wr2, wi;
00950                   const double *aa_ptr = aa.data () + jj * nn + jj;
00951                   const double *bb_ptr = bb.data () + jj * nn + jj;
00952                   F77_XFCN (dlag2, DLAG2,
00953                             (aa_ptr, nn, bb_ptr, nn, safmin,
00954                              scale1, scale2, wr1, wr2, wi));
00955 
00956 #ifdef DEBUG_EIG
00957                   std::cout << "dlag2 returns: scale1=" << scale1
00958                        << "\tscale2=" << scale2 << std::endl
00959                        << "\twr1=" << wr1 << "\twr2=" << wr2
00960                        << "\twi=" << wi << std::endl;
00961 #endif
00962 
00963                   // Just to be safe, check if it's a real pair.
00964                   if (wi == 0)
00965                     {
00966                       alphar(jj) = wr1;
00967                       alphai(jj) = 0;
00968                       betar(jj) = scale1;
00969                       alphar(jj+1) = wr2;
00970                       alphai(jj+1) = 0;
00971                       betar(jj+1) = scale2;
00972                     }
00973                   else
00974                     {
00975                       alphar(jj) = alphar(jj+1) = wr1;
00976                       alphai(jj) = -(alphai(jj+1) = wi);
00977                       betar(jj)  = betar(jj+1) = scale1;
00978                     }
00979                 }
00980 
00981               // Advance past this block.
00982               jj += zcnt;
00983             }
00984 
00985 #ifdef DEBUG_SORT
00986           std::cout << "qz: back from dsubsp: aa=" << std::endl;
00987           octave_print_internal (std::cout, aa, 0);
00988           std::cout << std::endl << "bb="  << std::endl;
00989           octave_print_internal (std::cout, bb, 0);
00990 
00991           if (compz == 'V')
00992             {
00993               std::cout << std::endl << "ZZ="  << std::endl;
00994               octave_print_internal (std::cout, ZZ, 0);
00995             }
00996           std::cout << std::endl << "qz: ndim=" << ndim << std::endl
00997                << "fail=" << fail << std::endl;
00998           std::cout << "alphar = " << std::endl;
00999           octave_print_internal (std::cout, (Matrix) alphar, 0);
01000           std::cout << std::endl << "alphai = " << std::endl;
01001           octave_print_internal (std::cout, (Matrix) alphai, 0);
01002           std::cout << std::endl << "beta = " << std::endl;
01003           octave_print_internal (std::cout, (Matrix) betar, 0);
01004           std::cout << std::endl;
01005 #endif
01006         }
01007     }
01008 
01009   // Compute generalized eigenvalues?
01010   ComplexColumnVector gev;
01011 
01012   if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4))
01013     {
01014       if (complex_case)
01015         {
01016           int cnt = 0;
01017 
01018           for (int ii = 0; ii < nn; ii++)
01019             cnt++;
01020 
01021           ComplexColumnVector tmp (cnt);
01022 
01023           cnt = 0;
01024           for (int ii = 0; ii < nn; ii++)
01025             tmp(cnt++) = xalpha(ii) / xbeta(ii);
01026 
01027           gev = tmp;
01028         }
01029       else
01030         {
01031 #ifdef DEBUG
01032           std::cout << "qz: computing generalized eigenvalues" << std::endl;
01033 #endif
01034 
01035           // Return finite generalized eigenvalues.
01036           int cnt = 0;
01037 
01038           for (int ii = 0; ii < nn; ii++)
01039             if (betar(ii) != 0)
01040               cnt++;
01041 
01042           ComplexColumnVector tmp (cnt);
01043 
01044           cnt = 0;
01045           for (int ii = 0; ii < nn; ii++)
01046             if (betar(ii) != 0)
01047               tmp(cnt++) = Complex(alphar(ii), alphai(ii))/betar(ii);
01048 
01049           gev = tmp;
01050         }
01051     }
01052 
01053   // Right, left eigenvector matrices.
01054   if (nargout >= 5)
01055     {
01056       // Which side to compute?
01057       char side = (nargout == 5 ? 'R' : 'B');
01058       // Compute all of them and backtransform
01059       char howmny = 'B';
01060       // Dummy pointer; select is not used.
01061       octave_idx_type *select = 0;
01062 
01063       if (complex_case)
01064         {
01065           CVL = CQ;
01066           CVR = CZ;
01067           ComplexRowVector cwork2 (2 * nn);
01068           RowVector rwork2 (8 * nn);
01069           octave_idx_type m;
01070 
01071           F77_XFCN (ztgevc, ZTGEVC,
01072                     (F77_CONST_CHAR_ARG2 (&side, 1),
01073                      F77_CONST_CHAR_ARG2 (&howmny, 1),
01074                      select, nn, caa.fortran_vec (), nn, cbb.fortran_vec (),
01075                      nn, CVL.fortran_vec (), nn, CVR.fortran_vec (), nn, nn,
01076                      m, cwork2.fortran_vec (), rwork2.fortran_vec (), info
01077                      F77_CHAR_ARG_LEN (1)
01078                      F77_CHAR_ARG_LEN (1)));
01079         }
01080       else
01081         {
01082 #ifdef DEBUG
01083           std::cout << "qz: computing  generalized eigenvectors" << std::endl;
01084 #endif
01085 
01086           VL = QQ;
01087           VR = ZZ;
01088           octave_idx_type m;
01089 
01090           F77_XFCN (dtgevc, DTGEVC,
01091                     (F77_CONST_CHAR_ARG2 (&side, 1),
01092                      F77_CONST_CHAR_ARG2 (&howmny, 1),
01093                      select, nn, aa.fortran_vec (), nn, bb.fortran_vec (),
01094                      nn, VL.fortran_vec (), nn, VR.fortran_vec (), nn, nn,
01095                      m, work.fortran_vec (), info
01096                      F77_CHAR_ARG_LEN (1)
01097                      F77_CHAR_ARG_LEN (1)));
01098 
01099           // Now construct the complex form of VV, WW.
01100           int jj = 0;
01101 
01102           while (jj < nn)
01103             {
01104               OCTAVE_QUIT;
01105 
01106               // See if real or complex eigenvalue.
01107 
01108               // Column increment; assume complex eigenvalue.
01109               int cinc = 2;
01110 
01111               if (jj == (nn-1))
01112                 // Single column.
01113                 cinc = 1;
01114               else if (aa(jj+1,jj) == 0)
01115                 cinc = 1;
01116 
01117               // Now copy the eigenvector (s) to CVR, CVL.
01118               if (cinc == 1)
01119                 {
01120                   for (int ii = 0; ii < nn; ii++)
01121                     CVR(ii,jj) = VR(ii,jj);
01122 
01123                   if (side == 'B')
01124                     for (int ii = 0; ii < nn; ii++)
01125                       CVL(ii,jj) = VL(ii,jj);
01126                 }
01127               else
01128                 {
01129                   // Double column; complex vector.
01130 
01131                   for (int ii = 0; ii < nn; ii++)
01132                     {
01133                       CVR(ii,jj) = Complex (VR(ii,jj), VR(ii,jj+1));
01134                       CVR(ii,jj+1) = Complex (VR(ii,jj), -VR(ii,jj+1));
01135                     }
01136 
01137                   if (side == 'B')
01138                     for (int ii = 0; ii < nn; ii++)
01139                       {
01140                         CVL(ii,jj) = Complex (VL(ii,jj), VL(ii,jj+1));
01141                         CVL(ii,jj+1) = Complex (VL(ii,jj), -VL(ii,jj+1));
01142                       }
01143                 }
01144 
01145               // Advance to next eigenvectors (if any).
01146               jj += cinc;
01147             }
01148         }
01149     }
01150 
01151   switch (nargout)
01152     {
01153     case 7:
01154       retval(6) = gev;
01155 
01156     case 6:
01157       // Return eigenvectors.
01158       retval(5) = CVL;
01159 
01160     case 5:
01161       // Return eigenvectors.
01162       retval(4) = CVR;
01163 
01164     case 4:
01165       if (nargin == 3)
01166         {
01167 #ifdef DEBUG
01168           std::cout << "qz: sort: retval(3) = gev = " << std::endl;
01169           octave_print_internal (std::cout, gev);
01170           std::cout << std::endl;
01171 #endif
01172           retval(3) = gev;
01173         }
01174       else
01175         {
01176           if (complex_case)
01177             retval(3) = CZ;
01178           else
01179             retval(3) = ZZ;
01180         }
01181 
01182     case 3:
01183       if (nargin == 3)
01184         retval(2) = CZ;
01185       else
01186         {
01187           if (complex_case)
01188             retval(2) = CQ.hermitian ();
01189           else
01190             retval(2) = QQ.transpose ();
01191         }
01192 
01193     case 2:
01194       {
01195         if (complex_case)
01196           {
01197 #ifdef DEBUG
01198             std::cout << "qz: retval (1) = cbb = " << std::endl;
01199             octave_print_internal (std::cout, cbb, 0);
01200             std::cout << std::endl << "qz: retval(0) = caa = " <<std::endl;
01201             octave_print_internal (std::cout, caa, 0);
01202             std::cout << std::endl;
01203 #endif
01204             retval(1) = cbb;
01205             retval(0) = caa;
01206           }
01207       else
01208         {
01209 #ifdef DEBUG
01210           std::cout << "qz: retval (1) = bb = " << std::endl;
01211           octave_print_internal (std::cout, bb, 0);
01212           std::cout << std::endl << "qz: retval(0) = aa = " <<std::endl;
01213           octave_print_internal (std::cout, aa, 0);
01214           std::cout << std::endl;
01215 #endif
01216           retval(1) = bb;
01217           retval(0) = aa;
01218         }
01219       }
01220       break;
01221 
01222 
01223     case 1:
01224     case 0:
01225 #ifdef DEBUG
01226       std::cout << "qz: retval(0) = gev = " << gev << std::endl;
01227 #endif
01228       retval(0) = gev;
01229       break;
01230 
01231     default:
01232       error ("qz: too many return arguments");
01233       break;
01234   }
01235 
01236 #ifdef DEBUG
01237   std::cout << "qz: exiting (at long last)" << std::endl;
01238 #endif
01239 
01240   return retval;
01241 }
01242 
01243 /*
01244 %!shared a, b, c
01245 %!  a = [1 2; 0 3];
01246 %!  b = [1 0; 0 0];
01247 %!  c = [0 1; 0 0];
01248 %!assert(qz (a,b), 1);
01249 %!assert(isempty (qz (a,c)));
01250 
01251 %% Exaple 7.7.3 in Golub & Van Loan
01252 %!test
01253 %! a = [ 10  1  2;
01254 %!        1  2 -1;
01255 %!        1  1  2];
01256 %! b = reshape(1:9,3,3);
01257 %! [aa, bb, q, z, v, w, lambda] = qz (a, b);
01258 %! sz = length(lambda);
01259 %! observed =  (b * v * diag ([lambda;0])) (:, 1:sz);
01260 %! assert ( (a*v) (:, 1:sz), observed, norm (observed) * 1e-14);
01261 %! observed = (diag ([lambda;0]) * w' * b) (1:sz, :);
01262 %! assert ( (w'*a) (1:sz, :) , observed, norm (observed) * 1e-13);
01263 %! assert (q * a * z, aa, norm (aa) * 1e-14);
01264 %! assert (q * b * z, bb, norm (bb) * 1e-14);
01265 
01266 %% FIXME: Still need a test for third form of calling qz
01267 
01268 */
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines