Mercurial > hg > octave-nkf
changeset 3185:9580887dd160
[project @ 1998-09-26 02:45:55 by jwe]
author | jwe |
---|---|
date | Sat, 26 Sep 1998 02:45:59 +0000 |
parents | 3988763ec9d3 |
children | edaa9a2d3d9c |
files | ChangeLog ROADMAP kpathsea/ChangeLog kpathsea/elt-dirs.c kpathsea/pathsearch.h liboctave/ChangeLog liboctave/Makefile.in liboctave/dMatrix.cc liboctave/dMatrix.h liboctave/mx-ext.h liboctave/oct-env.cc liboctave/pathsearch.cc liboctave/pathsearch.h src/ChangeLog src/DLD-FUNCTIONS/balance.cc src/DLD-FUNCTIONS/qz.cc src/Makefile.in src/defaults.cc src/help.cc src/load-save.cc |
diffstat | 20 files changed, 1065 insertions(+), 840 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog +++ b/ChangeLog @@ -1,6 +1,7 @@ Thu Sep 24 13:51:03 1998 John W. Eaton <jwe@bevo.che.wisc.edu> * configure.in (AC_OUTPUT): Add libcruft/ordered-qz to the list. + Delete libcruft/balgen and libcruft/eispack from the list. Mon Aug 31 12:07:02 1998 John W. Eaton <jwe@bevo.che.wisc.edu>
--- a/ROADMAP +++ b/ROADMAP @@ -14,16 +14,15 @@ kpathsea -- Karl Berry's path searching library libcruft -- various numerical libraries (mostly Fortran) - balgen * balancing subroutines for eigenvalue computations blas * basic linear algebra subroutines dassl * differential-algebraic system solver - eispack * subroutines for eigenvalue problems fftpack * subroutines for fast fourier transforms lapack * linear algebra package linpack * linear algebra package minpack * nonlinear equation solver misc * miscellaneous utilities odepack * odinary differential equation solver + ordered-qz * code for ordering eigenvalues for QZ factorization quadpack * subroutines for numerical integration ranlib * random number generators slatec-fn * various special function subroutines
--- a/kpathsea/ChangeLog +++ b/kpathsea/ChangeLog @@ -1,3 +1,8 @@ +Fri Sep 25 13:28:54 1998 John W. Eaton <jwe@bevo.che.wisc.edu> + + * elt-dirs.c (kpse_clear_dir_cache): New function. + * pathsearch.h (kpse_clear_dir_cache): Provide declaration. + Mon Aug 31 12:07:02 1998 John W. Eaton <jwe@bevo.che.wisc.edu> * config.sub: Accept armv4 everywhere arm is allowed.
--- a/kpathsea/elt-dirs.c +++ b/kpathsea/elt-dirs.c @@ -68,6 +68,33 @@ static cache_entry *the_cache = NULL; static unsigned cache_length = 0; +void +kpse_clear_dir_cache P1H(void) +{ + while (cache_length > 0) + { + str_llist_type elt = *the_cache[--cache_length].value; + + while (elt) + { + str_llist_type next = STR_LLIST_NEXT (*elt); + + string s = STR_LLIST (*elt); + + if (s) + free (s); + + free (elt); + + elt = next; + } + } + + if (the_cache) + free (the_cache); + + the_cache = NULL; +} /* Associate KEY with VALUE. We implement the cache as a simple linear list, since it's unlikely to ever be more than a dozen or so elements
--- a/kpathsea/pathsearch.h +++ b/kpathsea/pathsearch.h @@ -72,5 +72,8 @@ all the filenames (or NULL if none), instead of taking the first. */ extern string *kpse_all_path_search P2H(const_string path, const_string name); +/* Clear the directory cache. */ +extern void kpse_clear_dir_cache P1H(void); + #endif /* not KPATHSEA_PATHSEARCH_H */
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,5 +1,26 @@ +Fri Sep 25 14:26:44 1998 John W. Eaton <jwe@bevo.che.wisc.edu> + + * oct-env.cc (octave_env::do_get_home_directory): + If HOME can't be found, set it to "/". + (octave_env::do_get_user_name) + If user name can't be found, set it to "unknown". + (octave_env::do_get_host_name) + If host name can't be found, set it to "unknown". + + * pathsearch.h (dir_path::rehash): New function. + * pathsearch.cc (dir_path::init): Clear kpathsea's internal + diretcory cache before doing initialization. + Thu Sep 24 13:23:25 1998 John W. Eaton <jwe@bevo.che.wisc.edu> + * dMatrix.cc (Qzval): Delete. + (qzhes, qzit, qzval): Delete F77_FCN declarations. + * dMatrix.h (Qzval): Delete declaration. + + * dbleGEPBAL.h, dbleGEPBAL.cc: Delete. + * Makefile.in (MATRIX_INC, MATRIX_SRC): Delete them from the lists. + * mx-ext.h: Don't include dbleGEPBAL. + * lo-ieee.cc (octave_ieee_init): For now, use X_CAST instead of static_cast.
--- a/liboctave/Makefile.in +++ b/liboctave/Makefile.in @@ -31,7 +31,7 @@ CmplxCHOL.h CmplxDET.h CmplxHESS.h CmplxLU.h CmplxQR.h \ CmplxQRP.h CmplxSCHUR.h CmplxSVD.h EIG.h boolMatrix.h chMatrix.h \ dColVector.h dDiagMatrix.h dMatrix.h dRowVector.h dbleAEPBAL.h \ - dbleCHOL.h dbleDET.h dbleGEPBAL.h dbleHESS.h dbleLU.h dbleQR.h \ + dbleCHOL.h dbleDET.h dbleHESS.h dbleLU.h dbleQR.h \ dbleQRP.h dbleSCHUR.h dbleSVD.h MX_OP_INC := mx-cdm-cm.h mx-cdm-cs.h mx-cdm-dm.h mx-cdm-m.h \ @@ -65,9 +65,8 @@ CmplxHESS.cc CmplxLU.cc CmplxQR.cc CmplxQRP.cc CmplxSCHUR.cc \ CmplxSVD.cc EIG.cc MArray-misc.cc boolMatrix.cc chMatrix.cc \ dColVector.cc dDiagMatrix.cc dMatrix.cc dRowVector.cc \ - dbleAEPBAL.cc dbleCHOL.cc dbleDET.cc dbleGEPBAL.cc \ - dbleHESS.cc dbleLU.cc dbleQR.cc dbleQRP.cc dbleSCHUR.cc \ - dbleSVD.cc + dbleAEPBAL.cc dbleCHOL.cc dbleDET.cc dbleHESS.cc dbleLU.cc \ + dbleQR.cc dbleQRP.cc dbleSCHUR.cc dbleSVD.cc MX_OP_SRC := mx-cdm-cm.cc mx-cdm-cs.cc mx-cdm-dm.cc mx-cdm-m.cc \ mx-cdm-s.cc mx-cm-cdm.cc mx-cm-dm.cc mx-cm-m.cc mx-cm-s.cc \
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -99,17 +99,6 @@ int F77_FCN (xdlange, XDLANGE) (const char*, const int&, const int&, const double*, const int&, double*, double&); - - int F77_FCN (qzhes, QZHES) (const int&, const int&, double*, - double*, const long&, double*); - - int F77_FCN (qzit, QZIT) (const int&, const int&, double*, double*, - const double&, const long&, double*, - int&); - - int F77_FCN (qzval, QZVAL) (const int&, const int&, double*, - double*, double*, double*, double*, - const long&, double*); } // Matrix class. @@ -2935,94 +2924,6 @@ return retval; } -ComplexColumnVector -Qzval (const Matrix& a, const Matrix& b) -{ - ComplexColumnVector retval; - - int a_nr = a.rows(); - int a_nc = a.cols(); - - int b_nr = b.rows(); - int b_nc = b.cols(); - - if (a_nr == a_nc) - { - if (a_nr == b_nr && a_nc == b_nc) - { - if (a_nr != 0) - { - Matrix jnk (a_nr, a_nr, 0.0); - double *pjnk = jnk.fortran_vec (); - - ColumnVector alfr (a_nr); - double *palfr = alfr.fortran_vec (); - - ColumnVector alfi (a_nr); - double *palfi = alfi.fortran_vec (); - - ColumnVector beta (a_nr); - double *pbeta = beta.fortran_vec (); - - Matrix atmp = a; - double *pa = atmp.fortran_vec (); - - Matrix btmp = b; - double *pb = btmp.fortran_vec (); - - long matz = 0; - int info; - - // XXX FIXME ??? XXX - double eps = DBL_EPSILON; - - F77_FCN (qzhes, QZHES) (a_nr, a_nr, pa, pb, matz, pjnk); - - F77_FCN (qzit, QZIT) (a_nr, a_nr, pa, pb, eps, matz, pjnk, info); - - if (! info) - { - F77_FCN (qzval, QZVAL) (a_nr, a_nr, pa, pb, palfr, - palfi, pbeta, matz, pjnk); - - // Count and extract finite generalized eigenvalues. - - int cnt = 0; - - for (int i = 0; i < a_nr; i++) - if (beta(i) != 0) - cnt++; - - ComplexColumnVector cx (cnt); - - cnt = 0; - - for (int i = 0; i < a_nr; i++) - { - if (beta(i) != 0) - { - // Finite generalized eigenvalue. - - cx(cnt++) = Complex (alfr(i), alfi(i)) / beta(i); - } - } - - retval = cx; - } - else - (*current_liboctave_error_handler) - ("qzval: trouble in qzit, info = %d", info); - } - } - else - gripe_nonconformant ("qzval", a_nr, a_nc, b_nr, b_nc); - } - else - (*current_liboctave_error_handler) ("qzval: square matrices required"); - - return retval; -} - // matrix by matrix -> matrix operations Matrix
--- a/liboctave/dMatrix.h +++ b/liboctave/dMatrix.h @@ -234,8 +234,6 @@ extern Matrix Sylvester (const Matrix&, const Matrix&, const Matrix&); -extern ComplexColumnVector Qzval (const Matrix& a, const Matrix& b); - extern Matrix operator * (const Matrix& a, const Matrix& b); MS_CMP_OP_DECLS (Matrix, double)
--- a/liboctave/mx-ext.h +++ b/liboctave/mx-ext.h @@ -33,12 +33,6 @@ #include "dbleDET.h" #include "CmplxDET.h" -// Result of a GEP Balance operation -// Note: currenlty only do balancing on real data. Complex balancing -// done on magnitudes of complex data. - -#include "dbleGEPBAL.h" - // Result of a Cholesky Factorization #include "dbleCHOL.h"
--- a/liboctave/oct-env.cc +++ b/liboctave/oct-env.cc @@ -364,7 +364,7 @@ { octave_passwd pw = octave_passwd::getpwuid (octave_syscalls::getuid ()); - hd = pw ? pw.dir () : string ("I have no home!"); + hd = pw ? pw.dir () : string ("/"); } return hd; @@ -380,7 +380,7 @@ { octave_passwd pw = octave_passwd::getpwuid (octave_syscalls::getuid ()); - user_name = pw ? string ("I have no name!") : pw.name (); + user_name = pw ? string ("unknown") : pw.name (); } return user_name; @@ -398,7 +398,7 @@ int status = gethostname (hostname, 255); - host_name = (status < 0) ? "I have no host!" : hostname; + host_name = (status < 0) ? "unknown" : hostname; } return host_name;
--- a/liboctave/pathsearch.cc +++ b/liboctave/pathsearch.cc @@ -159,6 +159,8 @@ octave_kpathsea_initialized = true; } + kpse_clear_dir_cache (); + char *tmp = kpse_path_expand (p_orig.c_str ()); if (tmp) {
--- a/liboctave/pathsearch.h +++ b/liboctave/pathsearch.h @@ -71,6 +71,12 @@ static void set_program_name (const string&); + void rehash (void) + { + initialized = false; + init (); + } + private: // The colon separated list that we were given.
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,5 +1,28 @@ +Fri Sep 25 11:50:44 1998 John W. Eaton <jwe@bevo.che.wisc.edu> + + * load-save.cc (write_header): Rename from write_binary_header. + Also write header for Octave ASCII files. + + * load-save.cc (Fsave): Implement -append option. + + * defaults.cc (Frehash): New function. + +Fri Sep 25 11:50:44 1998 John W. Eaton <jwe@bevo.che.wisc.edu> + + * help.cc (help_from_info): Improve error message in case that + info doesn't work. + Thu Sep 24 10:48:12 1998 John W. Eaton <jwe@bevo.che.wisc.edu> + * Makefile.in (DLD_XSRC): Replace qzval.cc with qz.cc + + * DLD-FUNCTIONS/balance.cc: Update from A. S. Hodel + <scotte@eng.auburn.edu>. + + * DLD-FUNCTIONS/qz.cc: New file. + + * DLD-FUNCTIONS/qzval.cc: Delete. + * parse.y (plot_command1): Don't allow it to be empty. (plot_command): Handle simple `PLOT' and `PLOT ranges' as special cases here.
--- a/src/DLD-FUNCTIONS/balance.cc +++ b/src/DLD-FUNCTIONS/balance.cc @@ -42,24 +42,27 @@ extern "C" { - int F77_FCN( dggbal, DGGBAL) (const char* JOB, const int& N, - double* A, const int& LDA, double* B, const int& LDB, - int& ILO, int & IHI, double* LSCALE, - double* RSCALE, double* WORK, int& INFO, long ); + int F77_FCN (dggbal, DGGBAL) (const char* JOB, const int& N, + double* A, const int& LDA, double* B, + const int& LDB, int& ILO, int& IHI, + double* LSCALE, double* RSCALE, + double* WORK, int& INFO, long); - int F77_FCN( dggbak, DGGBAK) (const char* JOB, const char* SIDE, - const int& N, const int& ILO, const int& IHI, - double* LSCALE, double* RSCALE, int& M, - double* V, const int& LDV, int& INFO, long, long); + int F77_FCN (dggbak, DGGBAK) (const char* JOB, const char* SIDE, + const int& N, const int& ILO, + const int& IHI, double* LSCALE, + double* RSCALE, int& M, double* V, + const int& LDV, int& INFO, long, long); - int F77_FCN( zggbal, ZGGBAL) ( const char* JOB, const int& N, - Complex* A, const int& LDA, Complex* B, const int& LDB, - int& ILO, int & IHI, double* LSCALE, - double* RSCALE, double* WORK, int& INFO, long ); + int F77_FCN (zggbal, ZGGBAL) (const char* JOB, const int& N, + Complex* A, const int& LDA, Complex* B, + const int& LDB, int& ILO, int& IHI, + double* LSCALE, double* RSCALE, + double* WORK, int& INFO, long); } DEFUN_DLD (balance, args, nargout, - "AA = balance (A [, OPT]) or [[DD,] AA] = balance (A [, OPT])\n\ + "AA = balance (A [, OPT]) or [[DD,] AA] = balance (A [, OPT])\n\ \n\ generalized eigenvalue problem:\n\ \n\ @@ -77,7 +80,6 @@ \n\ [CC, DD, AA, BB] = balance (A, B, OPT) returns AA (BB) = CC*A*DD (CC*B*DD)") { - octave_value_list retval; int nargin = args.length (); @@ -89,7 +91,7 @@ } // determine if it's AEP or GEP - int AEPcase = (nargin == 1 ? 1 : args(1).is_string() ); + int AEPcase = nargin == 1 ? 1 : args(1).is_string (); string bal_job; // problem dimension @@ -99,181 +101,224 @@ if (arg_is_empty < 0) return retval; + if (arg_is_empty > 0) return octave_value_list (2, Matrix ()); if (nn != args(0).columns()) - { - gripe_square_matrix_required ("balance"); - return retval; - } + { + gripe_square_matrix_required ("balance"); + return retval; + } // Extract argument 1 parameter for both AEP and GEP. Matrix aa; ComplexMatrix caa; - if (args(0).is_complex_type ()) caa = args(0).complex_matrix_value (); - else aa = args(0).matrix_value (); - if (error_state) return retval; + + if (args(0).is_complex_type ()) + caa = args(0).complex_matrix_value (); + else + aa = args(0).matrix_value (); + + if (error_state) + return retval; // Treat AEP/GEP cases. - if(AEPcase) - { - // Algebraic eigenvalue problem. - if(nargin == 1) - bal_job = "B"; - else if(args(1).is_string()) - bal_job = args(1).string_value(); - // the next line should never execute, but better safe than sorry. - else error("balance: AEP argument 2 must be a string"); + if (AEPcase) + { + // Algebraic eigenvalue problem. - // balance the AEP - if (args(0).is_complex_type ()) - { - ComplexAEPBALANCE result (caa, bal_job); - - if (nargout == 0 || nargout == 1) - retval(0) = result.balanced_matrix (); + if (nargin == 1) + bal_job = "B"; + else if (args(1).is_string ()) + bal_job = args(1).string_value (); else - { - retval(1) = result.balanced_matrix (); - retval(0) = result.balancing_matrix (); - } - } - else - { - AEPBALANCE result (aa, bal_job); + { + error ("balance: AEP argument 2 must be a string"); + return retval; + } + + // balance the AEP + if (args(0).is_complex_type ()) + { + ComplexAEPBALANCE result (caa, bal_job); - if (nargout == 0 || nargout == 1) - retval(0) = result.balanced_matrix (); + if (nargout == 0 || nargout == 1) + retval(0) = result.balanced_matrix (); + else + { + retval(1) = result.balanced_matrix (); + retval(0) = result.balancing_matrix (); + } + } else - { - retval(1) = result.balanced_matrix (); - retval(0) = result.balancing_matrix (); - } + { + AEPBALANCE result (aa, bal_job); + + if (nargout == 0 || nargout == 1) + retval(0) = result.balanced_matrix (); + else + { + retval(1) = result.balanced_matrix (); + retval(0) = result.balancing_matrix (); + } + } } - } - // - // end of AEP case, now do GEP case else - { - // Generalized eigenvalue problem. - if(nargin == 2) - bal_job = "B"; - else if(args(2).is_string()) - bal_job = args(2).string_value(); - else error("balance: GEP argument 3 must be a string"); - - if( (nn != args(1).columns()) || (nn != args(1).rows() ) ) - { - gripe_nonconformant (); - return retval; - } - Matrix bb; - ComplexMatrix cbb; - if (args(1).is_complex_type ()) cbb = args(1).complex_matrix_value (); - else bb = args(1).matrix_value (); - if (error_state) return retval; - - // - // Both matrices loaded, now let's check what kind of arithmetic: - // first, declare variables used in both the real and complex case - int ilo, ihi, info; - RowVector lscale(nn), rscale(nn), work(6*nn); - char job = bal_job[0]; - static int complex_case = (args(0).is_complex_type() - || args(1).is_complex_type()); - - // now balance - if (complex_case) - { - if (args(0).is_real_type ()) caa = aa; - if (args(1).is_real_type ()) cbb = bb; - - F77_XFCN( zggbal, ZGGBAL, ( &job, nn, caa.fortran_vec(), nn, - cbb.fortran_vec(), nn, ilo, ihi, lscale.fortran_vec(), - rscale.fortran_vec(), work.fortran_vec(), info, 1L)); - } - else // real matrices case { - F77_XFCN( dggbal, DGGBAL, (&job, nn, aa.fortran_vec(), - nn, bb.fortran_vec() , nn, ilo, ihi, lscale.fortran_vec(), - rscale.fortran_vec(), work.fortran_vec(), info , 1L)); - - if(f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in balance GEP"); - } - - // Since we just want the balancing matrices, we can use dggbal - // for both the real and complex cases; - Matrix Pl(nn,nn), Pr(nn,nn); - for(int ii=0; ii < nn ; ii++) - for( int jj=0; jj < nn ; jj++) - Pl(ii,jj) = Pr(ii,jj) = (ii == jj ? 1.0 : 0.0); + // Generalized eigenvalue problem. + if (nargin == 2) + bal_job = "B"; + else if (args(2).is_string ()) + bal_job = args(2).string_value (); + else + { + error ("balance: GEP argument 3 must be a string"); + return retval; + } + + if ((nn != args(1).columns ()) || (nn != args(1).rows ())) + { + gripe_nonconformant (); + return retval; + } + + Matrix bb; + ComplexMatrix cbb; + + if (args(1).is_complex_type ()) + cbb = args(1).complex_matrix_value (); + else + bb = args(1).matrix_value (); + + if (error_state) + return retval; + + // Both matrices loaded, now let's check what kind of arithmetic: + // first, declare variables used in both the real and complex case + + int ilo, ihi, info; + RowVector lscale(nn), rscale(nn), work(6*nn); + char job = bal_job[0]; + + static int complex_case + = (args(0).is_complex_type () || args(1).is_complex_type ()); + + // now balance + if (complex_case) + { + if (args(0).is_real_type ()) + caa = aa; + + if (args(1).is_real_type ()) + cbb = bb; - // left first - F77_XFCN( dggbak, DGGBAK, (&job, "L", - nn, ilo, ihi, lscale.fortran_vec(), - rscale.fortran_vec(), nn, Pl.fortran_vec(), - nn, info, 1L, 1L)); - - if(f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in balance GEP(L)"); - - // then right - F77_XFCN(dggbak, DGGBAK, (&job, "R", - nn, ilo, ihi, lscale.fortran_vec(), - rscale.fortran_vec(), nn, Pr.fortran_vec(), - nn, info, 1L, 1L)); - if(f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in balance GEP(R)"); + F77_XFCN (zggbal, ZGGBAL, + (&job, nn, caa.fortran_vec(), nn, + cbb.fortran_vec(), nn, ilo, ihi, + lscale.fortran_vec(), rscale.fortran_vec(), + work.fortran_vec(), info, 1L)); + + if (f77_exception_encountered) + { + error ("unrecoverable error in balance GEP"); + return retval; + } + } + else + { + // real matrices case - switch (nargout) - { - case 0: - case 1: - warning ("balance: used GEP, should have two output arguments"); - if(complex_case) - retval(0) = caa; - else - retval(0) = aa; - break; + F77_XFCN (dggbal, DGGBAL, + (&job, nn, aa.fortran_vec(), nn, bb.fortran_vec(), + nn, ilo, ihi, lscale.fortran_vec(), + rscale.fortran_vec(), work.fortran_vec(), info, 1L)); + + if (f77_exception_encountered) + { + error ("unrecoverable error in balance GEP"); + return retval; + } + } + + // Since we just want the balancing matrices, we can use dggbal + // for both the real and complex cases. + + Matrix Pl(nn,nn), Pr(nn,nn); + + for (int ii = 0; ii < nn; ii++) + for (int jj = 0; jj < nn; jj++) + Pl(ii,jj) = Pr(ii,jj) = (ii == jj ? 1.0 : 0.0); + + // left first + F77_XFCN (dggbak, DGGBAK, + (&job, "L", nn, ilo, ihi, lscale.fortran_vec(), + rscale.fortran_vec(), nn, Pl.fortran_vec(), + nn, info, 1L, 1L)); + + if (f77_exception_encountered) + { + error ("unrecoverable error in balance GEP(L)"); + return retval; + } + + // then right + F77_XFCN (dggbak, DGGBAK, + (&job, "R", nn, ilo, ihi, lscale.fortran_vec(), + rscale.fortran_vec(), nn, Pr.fortran_vec(), + nn, info, 1L, 1L)); + + if (f77_exception_encountered) + { + error ("unrecoverable error in balance GEP(R)"); + return retval; + } - case 2: - if(complex_case) - { - retval(1) = cbb; - retval(0) = caa; - } - else - { - retval(1) = bb; - retval(0) = aa; - } - break; + switch (nargout) + { + case 0: + case 1: + warning ("balance: used GEP, should have two output arguments"); + if (complex_case) + retval(0) = caa; + else + retval(0) = aa; + break; - case 4: - if(complex_case) - { - retval(3) = cbb; - retval(2) = caa; - } - else - { - retval(3) = bb; - retval(2) = aa; - } - retval(1) = Pr; - retval(0) = Pl.transpose(); // so that aa_bal = cc*aa*dd, etc. - break; + case 2: + if (complex_case) + { + retval(1) = cbb; + retval(0) = caa; + } + else + { + retval(1) = bb; + retval(0) = aa; + } + break; - default: - error ("balance: invalid number of output arguments"); - break; + case 4: + if (complex_case) + { + retval(3) = cbb; + retval(2) = caa; + } + else + { + retval(3) = bb; + retval(2) = aa; + } + retval(1) = Pr; + retval(0) = Pl.transpose (); // so that aa_bal = cc*aa*dd, etc. + break; + + default: + error ("balance: invalid number of output arguments"); + break; + } } - } + return retval; }
--- a/src/DLD-FUNCTIONS/qz.cc +++ b/src/DLD-FUNCTIONS/qz.cc @@ -47,61 +47,77 @@ #include "oct-map.h" #include "ov.h" #include "pager.h" -#if defined(DEBUG) || defined(DEBUG_SORT) +#if defined (DEBUG) || defined (DEBUG_SORT) #include "pr-output.h" #endif #include "symtab.h" #include "utils.h" #include "variables.h" -typedef int (*sort_function) (const int& LSIZE, const double& ALPHA, - const double& BETA, const double& S, const double& P); +typedef int (*sort_function) (const int& LSIZE, const double& ALPHA, + const double& BETA, const double& S, + const double& P); extern "C" { - int F77_FCN( dggbal, DGGBAL) (const char* JOB, const int& N, - double* A, const int& LDA, double* B, const int& LDB, - int& ILO, int & IHI, double* LSCALE, - double* RSCALE, double* WORK, int& INFO, long ); + int F77_FCN (dggbal, DGGBAL) (const char* JOB, const int& N, + double* A, const int& LDA, double* B, + const int& LDB, int& ILO, int& IHI, + double* LSCALE, double* RSCALE, + double* WORK, int& INFO, long); - int F77_FCN( dggbak, DGGBAK) (const char* JOB, const char* SIDE, - const int& N, const int& ILO, const int& IHI, - double* LSCALE, double* RSCALE, int& M, - double* V, const int& LDV, int& INFO, long, long); + int F77_FCN (dggbak, DGGBAK) (const char* JOB, const char* SIDE, + const int& N, const int& ILO, + const int& IHI, double* LSCALE, + double* RSCALE, int& M, double* V, + const int& LDV, int& INFO, long, long); - int F77_FCN( dgghrd, DGGHRD) ( const char* COMPQ, const char* COMPZ, - const int& N, const int& ILO, const int& IHI, double* A, const int& LDA, - double* B, const int& LDB, double* Q, const int& LDQ, double* Z, - const int& LDZ, int& INFO, const long, const long); + int F77_FCN (dgghrd, DGGHRD) (const char* COMPQ, const char* COMPZ, + const int& N, const int& ILO, + const int& IHI, double* A, + const int& LDA, double* B, + const int& LDB, double* Q, + const int& LDQ, double* Z, + const int& LDZ, int& INFO, long, long); - int F77_FCN( dhgeqz, DHGEQZ) ( const char* JOB, const char* COMPQ, - const char* COMPZ, const int& N, const int& ILO, const int& IHI, - double* A, const int& LDA, double* B, const int& LDB, - double* ALPHAR, double* ALPHAI, double* BETA, double* Q, - const int& LDQ, double* Z, const int& LDZ, double* WORK, - const int& LWORK, int& INFO, const long, const long, const long ); + int F77_FCN (dhgeqz, DHGEQZ) (const char* JOB, const char* COMPQ, + const char* COMPZ, const int& N, + const int& ILO, const int& IHI, + double* A, const int& LDA, double* B, + const int& LDB, double* ALPHAR, + double* ALPHAI, double* BETA, double* Q, + const int& LDQ, double* Z, + const int& LDZ, double* WORK, + const int& LWORK, int& INFO, + long, long, long); - int F77_FCN( dlag2, DLAG2) ( double* A, const int& LDA, double* B, - const int& LDB, const double& SAFMIN, double& SCALE1, double& SCALE2, - double& WR1, double& WR2, double& WI ); + int F77_FCN (dlag2, DLAG2) (double* A, const int& LDA, double* B, + const int& LDB, const double& SAFMIN, + double& SCALE1, double& SCALE2, + double& WR1, double& WR2, double& WI); // Van Dooren's code (netlib.org: toms/590) for reordering // GEP. Only processes Z, not Q. - int F77_FCN( dsubsp, DSUBSP) ( const int& NMAX, const int& N, double* A, - double* B, double* Z, sort_function, const double& EPS, - int& NDIM, int& FAIL, int* IND); + int F77_FCN (dsubsp, DSUBSP) (const int& NMAX, const int& N, double* A, + double* B, double* Z, sort_function, + const double& EPS, int& NDIM, int& FAIL, + int* IND); // documentation for DTGEVC incorrectly states that VR, VL are // complex*16; they are declared in DTGEVC as double precision // (probably a cut and paste problem fro ZTGEVC) - int F77_FCN( dtgevc, DTGEVC) ( const char* SIDE, const char* HOWMNY, - int* SELECT, const int& N, double* A, const int& LDA, double* B, - const int& LDB, double* VL, const int& LDVL, double* VR, - const int& LDVR, const int& MM, int& M, double* WORK, int& INFO, - long, long ); + int F77_FCN (dtgevc, DTGEVC) (const char* SIDE, const char* HOWMNY, + int* SELECT, const int& N, double* A, + const int& LDA, double* B, + const int& LDB, double* VL, + const int& LDVL, double* VR, + const int& LDVR, const int& MM, + int& M, double* WORK, int& INFO, + long, long); - int F77_FCN ( xdlamch, XDLAMCH) (const char* cmach, double& retval, long); - int F77_FCN ( xdlange, XDLANGE) (const char*, const int&, + int F77_FCN (xdlamch, XDLAMCH) (const char* cmach, double& retval, long); + + int F77_FCN (xdlange, XDLANGE) (const char*, const int&, const int&, const double*, const int&, double*, double&); } @@ -114,64 +130,79 @@ // fcrhp: real(lambda) >= 0 // folhp: real(lambda) < 0 -static int fcrhp(const int& lsize, const double& alpha, - const double& beta, const double& s, const double& p) +static int +fcrhp (const int& lsize, const double& alpha, + const double& beta, const double& s, const double&) { - if(lsize == 1) + if (lsize == 1) return (alpha*beta >= 0 ? 1 : -1); - else + else return (s >= 0 ? 1 : -1); } -static int fin(const int& lsize, const double& alpha, - const double& beta, const double& s, const double& p) + +static int +fin (const int& lsize, const double& alpha, + const double& beta, const double&, const double& p) { int retval; - if(lsize == 1) - retval = (fabs(alpha) < fabs(beta) ? 1 : -1); - else - retval = (fabs(p) < 1 ? 1 : -1); - #ifdef DEBUG + if (lsize == 1) + retval = (fabs (alpha) < fabs (beta) ? 1 : -1); + else + retval = (fabs (p) < 1 ? 1 : -1); + +#ifdef DEBUG cout << "qz: fin: retval=" << retval << endl; - #endif +#endif + return retval; } -static int folhp(const int& lsize, const double& alpha, - const double& beta, const double& s, const double& p) + +static int +folhp (const int& lsize, const double& alpha, + const double& beta, const double& s, const double&) { - if(lsize == 1) + if (lsize == 1) return (alpha*beta < 0 ? 1 : -1); - else + else return (s < 0 ? 1 : -1); } -static int fout(const int& lsize, const double& alpha, - const double& beta, const double& s, const double& p) + +static int +fout (const int& lsize, const double& alpha, + const double& beta, const double&, const double& p) { - if(lsize == 1) - return (fabs(alpha) >= fabs(beta) ? 1 : -1); - else - return (fabs(p) >= 1 ? 1 : -1); + if (lsize == 1) + return (fabs (alpha) >= fabs (beta) ? 1 : -1); + else + return (fabs (p) >= 1 ? 1 : -1); } DEFUN_DLD (qz, args, nargout, -"Usage: lambda = qz(A,B) form [1]\n\ - [AA,BB,Q,Z{,V,W,lambda}] = qz(A,B) form [2]\n\ - [AA,BB,Z{,lambda}] = qz(A,B,opt) form [3]\n\ +"Usage:\n\ + + lambda = qz (A, B) form [1]\n\ + [AA, BB, Q, Z {, V, W, lambda}] = qz (A, B) form [2]\n\ + [AA, BB, Z{, lambda}] = qz (A, B, opt) form [3]\n\ +\n\ Generalized eigenvalue problem A x = s B x \n\ - +\n\ Form [1]: Computes the generalized eigenvalues lambda of (A - sB).\n\ +\n\ Form [2]: Computes qz decomposition, generalized eigenvectors, and \n\ - generalized eigenvalues of (A - sB)\n\ - A V = B V diag(lambda)\n\ - W' A = diag(lambda) W' B\n\ - AA = Q'*A*Z, BB = Q'*B*Z with Q, Z orthogonal (unitary)= I\n\ + generalized eigenvalues of (A - sB)\n\ + A V = B V diag (lambda)\n\ + W' A = diag (lambda) W' B\n\ + AA = Q'*A*Z, BB = Q'*B*Z with Q, Z orthogonal (unitary)= I\n\ +\n\ Form [3]: As in form [2], but allows ordering of generalized eigenpairs\n\ - for (e.g.) solution of discrete time algebraic Riccati equations.\n\ - Form 3 is not available for complex matrices and does not compute\n\ - the generalized eigenvectors V, W, nor the orthogonal matrix Q.\n\ + for (e.g.) solution of discrete time algebraic Riccati equations.\n\ + Form 3 is not available for complex matrices and does not compute\n\ + the generalized eigenvectors V, W, nor the orthogonal matrix Q.\n\ \n\ - opt: for ordering eigenvalues of the GEP pencil. The leading block\n\ - of the revised pencil contains all eigenvalues that satisfy:\n\ + opt: for ordering eigenvalues of the GEP pencil. The leading block\n\ + of the revised pencil contains all eigenvalues that satisfy:\n\ +\n\ \"N\" = unordered (default) \n\ \"S\" = small: leading block has all |lambda| <=1 \n\ \"B\" = big: leading block has all |lambda >= 1 \n\ @@ -183,586 +214,705 @@ Note: Permutation balancing is performed, but not scaling (see balance)\n\ Order of output arguments was selected for compatibility with MATLAB\n\ \n\ -See also: balance, dare, eig, schur\n") +See also: balance, dare, eig, schur") { octave_value_list retval; int nargin = args.length (); - #ifdef DEBUG +#ifdef DEBUG cout << "qz: nargin = " << nargin << ", nargout = " << nargout << endl; - #endif +#endif - if (nargin < 2 || nargin > 3 || nargout > 7 ) - { - print_usage ("qz"); - return retval; - } - else if(nargin == 3 && (nargout < 3 || nargout > 4)) - { - error("qz: Illegal number of output arguments for form [3] call"); - } + if (nargin < 2 || nargin > 3 || nargout > 7) + { + print_usage ("qz"); + return retval; + } + else if (nargin == 3 && (nargout < 3 || nargout > 4)) + { + error ("qz: Illegal number of output arguments for form [3] call"); + return retval; + } - #ifdef DEBUG +#ifdef DEBUG cout << "qz: determine ordering option" << endl; - #endif +#endif // Determine ordering option string ord_job; static double safmin; - if(nargin == 2) + + if (nargin == 2) ord_job = "N"; - else if( !args(2).is_string() ) - error("qz: argument 3 must be a string"); - else - { - ord_job = args(2).string_value(); - if(ord_job[0] != 'N' && ord_job[0] != 'S' && ord_job[0] != 'B' - && ord_job[0] != '+' && ord_job[0] != '-') - error("qz: illegal order option"); + else if (!args(2).is_string ()) + { + error ("qz: argument 3 must be a string"); + return retval; + } + else + { + ord_job = args(2).string_value (); - // overflow constant required by dlag2 - F77_XFCN ( xdlamch, XDLAMCH, ("S", safmin, 1L)); - - #ifdef DEBUG_EIG - cout << "qz: initial value of safmin=" << setiosflags(ios::scientific) - << safmin << endl; - #endif + if (ord_job[0] != 'N' + && ord_job[0] != 'S' + && ord_job[0] != 'B' + && ord_job[0] != '+' + && ord_job[0] != '-') + { + error ("qz: illegal order option"); + return retval; + } + + // overflow constant required by dlag2 + F77_FCN (xdlamch, XDLAMCH) ("S", safmin, 1L); - // some machines (e.g., DEC alpha) get safmin = 0; - // for these, use eps instead to avoid problems in dlag2 - if(safmin == 0) - { - #ifdef DEBUG_EIG - cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << endl; - #endif +#ifdef DEBUG_EIG + cout << "qz: initial value of safmin=" << setiosflags (ios::scientific) + << safmin << endl; +#endif - F77_XFCN ( xdlamch, XDLAMCH, ("E", safmin, 1L)); + // some machines (e.g., DEC alpha) get safmin = 0; + // for these, use eps instead to avoid problems in dlag2 + if (safmin == 0) + { +#ifdef DEBUG_EIG + cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << endl; +#endif - #ifdef DEBUG_EIG - cout << "qz: safmin set to " << setiosflags(ios::scientific) - << safmin << endl; - #endif + F77_FCN (xdlamch, XDLAMCH) ("E", safmin, 1L); + +#ifdef DEBUG_EIG + cout << "qz: safmin set to " << setiosflags (ios::scientific) + << safmin << endl; +#endif + } } - } - #ifdef DEBUG +#ifdef DEBUG cout << "qz: check argument 1" << endl; - #endif +#endif // Argument 1: check if it's o.k. dimensioned - int nn = args(0).rows(); - - #ifdef DEBUG - cout << "argument 1 dimensions: (" << nn << "," << args(0).columns() << ")" - << endl; - #endif - int arg_is_empty = empty_arg ("qz", nn, args(0).columns()); + int nn = args(0).rows (); + +#ifdef DEBUG + cout << "argument 1 dimensions: (" << nn << "," << args(0).columns () << ")" + << endl; +#endif + + int arg_is_empty = empty_arg ("qz", nn, args(0).columns ()); + if (arg_is_empty < 0) - { - gripe_empty_arg("qz: parameter 1",0); - return retval; - } + { + gripe_empty_arg ("qz: parameter 1", 0); + return retval; + } else if (arg_is_empty > 0) - { - gripe_empty_arg("qz: parameter 1; continuing",0); - return octave_value_list (2, Matrix ()); - } - else if (args(0).columns() != nn) - { - gripe_square_matrix_required ("qz"); - return retval; - } + { + gripe_empty_arg ("qz: parameter 1; continuing", 0); + return octave_value_list (2, Matrix ()); + } + else if (args(0).columns () != nn) + { + gripe_square_matrix_required ("qz"); + return retval; + } // Argument 1: dimensions look good; get the value Matrix aa; ComplexMatrix caa; - if (args(0).is_complex_type ()) + + if (args(0).is_complex_type ()) caa = args(0).complex_matrix_value (); - else + else aa = args(0).matrix_value (); - if (error_state) + + if (error_state) return retval; - #ifdef DEBUG +#ifdef DEBUG cout << "qz: check argument 2" << endl; - #endif +#endif // Extract argument 2 (bb, or cbb if complex) - if( (nn != args(1).columns()) || (nn != args(1).rows() )) - { - gripe_nonconformant (); - return retval; - } + if ((nn != args(1).columns ()) || (nn != args(1).rows ())) + { + gripe_nonconformant (); + return retval; + } + Matrix bb; ComplexMatrix cbb; - if (args(1).is_complex_type ()) + + if (args(1).is_complex_type ()) cbb = args(1).complex_matrix_value (); else bb = args(1).matrix_value (); - if (error_state) + + if (error_state) return retval; // Both matrices loaded, now let's check what kind of arithmetic: //declared static to avoid compiler warnings about long jumps, vforks. - static int complex_case - = (args(0).is_complex_type() || args(1).is_complex_type()); + + static int complex_case + = (args(0).is_complex_type () || args(1).is_complex_type ()); - if(nargin == 3 && complex_case) - error("qz: cannot re-order complex qz decomposition."); + if (nargin == 3 && complex_case) + { + error ("qz: cannot re-order complex qz decomposition."); + return retval; + } // first, declare variables used in both the real and complex case Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn); RowVector alphar(nn), alphai(nn), betar(nn); - ComplexMatrix CQ(nn,nn), CZ(nn,nn),CVR(nn,nn),CVL(nn,nn); + ComplexMatrix CQ(nn,nn), CZ(nn,nn), CVR(nn,nn), CVL(nn,nn); int ilo, ihi, info; - char compq = (nargout >= 3 ? 'V' : 'N'), - compz = (nargout >= 4 ? 'V' : 'N'); + char compq = (nargout >= 3 ? 'V' : 'N'); + char compz = (nargout >= 4 ? 'V' : 'N'); - // initialize Q,Z to identity if we need either of them - if(compq == 'V' || compz == 'V') - for(int ii=0; ii < nn ; ii++) - for( int jj=0; jj < nn ; jj++) + // initialize Q, Z to identity if we need either of them + if (compq == 'V' || compz == 'V') + for (int ii = 0; ii < nn; ii++) + for (int jj = 0; jj < nn; jj++) QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0); - // always perform permutation balancing + // always perform permutation balancing char bal_job = 'P'; RowVector lscale(nn), rscale(nn), work(6*nn); - if(complex_case) - error("Complex case not implemented yet"); + if (complex_case) + { + error ("Complex case not implemented yet"); + return retval; + } else - { - #ifdef DEBUG - if(compq == 'V') - cout << "qz: performing balancing; QQ=" << endl << QQ << endl; - #endif + { +#ifdef DEBUG + if (compq == 'V') + cout << "qz: performing balancing; QQ=" << endl << QQ << endl; +#endif - F77_XFCN( dggbal, DGGBAL, (&bal_job, nn, aa.fortran_vec(), - nn, bb.fortran_vec() , nn, ilo, ihi, lscale.fortran_vec(), - rscale.fortran_vec(), work.fortran_vec(), info , 1L)); - if(f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in qz(bal)"); - } + F77_XFCN (dggbal, DGGBAL, + (&bal_job, nn, aa.fortran_vec(), nn, bb.fortran_vec(), + nn, ilo, ihi, lscale.fortran_vec(), + rscale.fortran_vec(), work.fortran_vec(), info, 1L)); + + if (f77_exception_encountered) + { + error ("unrecoverable error in qz (bal)"); + return retval; + } + } // Since we just want the balancing matrices, we can use dggbal // for both the real and complex cases; // left first - if(compq == 'V') - { - F77_XFCN( dggbak, DGGBAK, (&bal_job, "L", - nn, ilo, ihi, lscale.fortran_vec(), - rscale.fortran_vec(), nn, QQ.fortran_vec(), - nn, info, 1L, 1L)); + + if (compq == 'V') + { + F77_XFCN (dggbak, DGGBAK, + (&bal_job, "L", nn, ilo, ihi, lscale.fortran_vec(), + rscale.fortran_vec(), nn, QQ.fortran_vec(), + nn, info, 1L, 1L)); - #ifdef DEBUG - if(compq == 'V') cout << "qz: balancing done; QQ=" << endl << QQ << endl; - #endif +#ifdef DEBUG + if (compq == 'V') + cout << "qz: balancing done; QQ=" << endl << QQ << endl; +#endif - if(f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in qz(bal-L)"); + if (f77_exception_encountered) + { + error ("unrecoverable error in qz (bal-L)"); + return retval; + } } // then right - if(compz == 'V') - { - F77_XFCN(dggbak, DGGBAK, (&bal_job, "R", - nn, ilo, ihi, lscale.fortran_vec(), - rscale.fortran_vec(), nn, ZZ.fortran_vec(), - nn, info, 1L, 1L)); + if (compz == 'V') + { + F77_XFCN (dggbak, DGGBAK, (&bal_job, "R", + nn, ilo, ihi, lscale.fortran_vec(), + rscale.fortran_vec(), nn, ZZ.fortran_vec(), + nn, info, 1L, 1L)); - #ifdef DEBUG - if(compz == 'V') cout << "qz: balancing done; ZZ=" << endl << ZZ << endl; - #endif +#ifdef DEBUG + if (compz == 'V') + cout << "qz: balancing done; ZZ=" << endl << ZZ << endl; +#endif - if(f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in qz(bal-R)"); - } + if (f77_exception_encountered) + { + error ("unrecoverable error in qz (bal-R)"); + return retval; + } + } static char qz_job; qz_job = (nargout < 2 ? 'E' : 'S'); + if (complex_case) - { - // complex case - if (args(0).is_real_type ()) caa = aa; - if (args(1).is_real_type ()) cbb = bb; - if(compq == 'V') CQ = QQ; - if(compz == 'V') CZ = ZZ; - error("complex case not done yet"); - } + { + // complex case + if (args(0).is_real_type ()) + caa = aa; + + if (args(1).is_real_type ()) + cbb = bb; + + if (compq == 'V') + CQ = QQ; + + if (compz == 'V') + CZ = ZZ; + + error ("complex case not done yet"); + return retval; + } else // real matrices case - { - #ifdef DEBUG - cout << "qz: peforming qr decomposition of bb" << endl; - #endif + { +#ifdef DEBUG + cout << "qz: peforming qr decomposition of bb" << endl; +#endif - // compute the QR factorization of bb - QR bqr(bb); + // compute the QR factorization of bb + QR bqr (bb); + +#ifdef DEBUG + cout << "qz: qr (bb) done; now peforming qz decomposition" << endl; +#endif - #ifdef DEBUG - cout << "qz: qr(bb) done; now peforming qz decomposition" << endl; - #endif + bb = bqr.R (); + +#ifdef DEBUG + cout << "qz: extracted bb" << endl; +#endif - bb = bqr.R(); - #ifdef DEBUG - cout << "qz: extracted bb" << endl; - #endif + aa = (bqr.Q ()).transpose ()*aa; + +#ifdef DEBUG + cout << "qz: updated aa " << endl; + cout << "bqr.Q () = " << endl << bqr.Q () << endl; - aa = (bqr.Q()).transpose()*aa; - #ifdef DEBUG - cout << "qz: updated aa " << endl; - cout << "bqr.Q () = " << endl << bqr.Q () << endl; - if(compq == 'V') cout << "QQ =" << QQ << endl; - #endif + if (compq == 'V') + cout << "QQ =" << QQ << endl; +#endif - if(compq == 'V') QQ = QQ*bqr.Q(); + if (compq == 'V') + QQ = QQ*bqr.Q (); - #ifdef DEBUG - cout << "qz: precursors done..." << endl; - #endif +#ifdef DEBUG + cout << "qz: precursors done..." << endl; +#endif +#ifdef DEBUG + cout << "qz: compq = " << compq << ", compz = " << compz << endl; +#endif - #ifdef DEBUG - cout << "qz: compq = " << compq << ", compz = " << compz << endl; - #endif + // reduce to generalized hessenberg form + F77_XFCN (dgghrd, DGGHRD, + (&compq, &compz, nn, ilo, ihi, aa.fortran_vec(), + nn, bb.fortran_vec(), nn, QQ.fortran_vec(), nn, + ZZ.fortran_vec(), nn, info, 1L, 1L)); - // reduce to generalized hessenberg form - F77_XFCN( dgghrd, DGGHRD, (&compq, &compz, nn, ilo, ihi, aa.fortran_vec(), - nn, bb.fortran_vec(), nn, QQ.fortran_vec(), nn, ZZ.fortran_vec(), - nn, info,1L,1L)); - if(f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in qz(dgghrd)"); + if (f77_exception_encountered) + { + error ("unrecoverable error in qz (dgghrd)"); + return retval; + } - // check if just computing generalized eigenvalues or if we're - // actually computing the decomposition + // check if just computing generalized eigenvalues or if we're + // actually computing the decomposition - // reduce to generalized Schur form - F77_XFCN( dhgeqz, DHGEQZ, ( &qz_job, &compq, &compz, nn, ilo, ihi, - aa.fortran_vec(), nn, bb.fortran_vec(), nn, alphar.fortran_vec(), - alphai.fortran_vec(), betar.fortran_vec(), QQ.fortran_vec(), - nn, ZZ.fortran_vec(), nn, work.fortran_vec(), nn, info, 1L, 1L, 1L)); - if(f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in qz(dhgeqz)"); + // reduce to generalized Schur form + F77_XFCN (dhgeqz, DHGEQZ, + (&qz_job, &compq, &compz, nn, ilo, ihi, + aa.fortran_vec(), nn, bb.fortran_vec(), nn, + alphar.fortran_vec(), alphai.fortran_vec(), + betar.fortran_vec(), QQ.fortran_vec(), nn, + ZZ.fortran_vec(), nn, work.fortran_vec(), nn, info, + 1L, 1L, 1L)); - } + if (f77_exception_encountered) + { + error ("unrecoverable error in qz (dhgeqz)"); + return retval; + } + } // order the QZ decomposition? - if(ord_job[0] != 'N') - { - if(complex_case) // probably not needed, but better be safe - error("qz: cannot re-order complex qz decomposition."); - - else + if (ord_job[0] != 'N') { - #ifdef DEBUG_SORT - cout << "qz: ordering eigenvalues: ord_job = " << ord_job[0] << endl; - #endif + if (complex_case) + { + // probably not needed, but better be safe + error ("qz: cannot re-order complex qz decomposition."); + return retval; + } + else + { +#ifdef DEBUG_SORT + cout << "qz: ordering eigenvalues: ord_job = " << ord_job[0] << endl; +#endif - // declared static to avoid vfork/long jump compiler complaints - static sort_function sort_test; - sort_test = NULL; - - switch(ord_job[0]) - { - case 'S': - sort_test = &fin; - break; - case 'B': - sort_test = &fout; - break; - case '+': - sort_test = &fcrhp; - break; - case '-': - sort_test = &folhp; - break; - default: // this should never happen - error("qz: illegal order option"); - } - - int ndim, fail, ind[nn]; - double inf_norm; - F77_XFCN (xdlange, XDLANGE, ("I", nn, nn, aa.fortran_vec (), nn, - work.fortran_vec (), inf_norm)); - - double eps = DBL_EPSILON*inf_norm*nn; - - #ifdef DEBUG_SORT - cout << "qz: calling dsubsp: aa=" << endl; - octave_print_internal(cout,aa,0); - cout << endl << "bb=" << endl; - octave_print_internal(cout,bb,0); - if(compz == 'V') - { - cout << endl << "ZZ=" << endl; - octave_print_internal(cout,ZZ,0); - } - cout << endl; - cout << "alphar = " << endl; - octave_print_internal(cout,(Matrix) alphar,0); - cout << endl << "alphai = " << endl; - octave_print_internal(cout,(Matrix) alphai,0); - cout << endl << "beta = " << endl; - octave_print_internal(cout,(Matrix) betar,0); - cout << endl; - #endif - - F77_XFCN( dsubsp, DSUBSP, (nn,nn,aa.fortran_vec(), bb.fortran_vec(), - ZZ.fortran_vec(), sort_test, eps, ndim, fail, ind)); + // declared static to avoid vfork/long jump compiler complaints + static sort_function sort_test; + sort_test = NULL; - #ifdef DEBUG - cout << "qz: back from dsubsp: aa=" << endl; - octave_print_internal(cout,aa,0); - cout << endl << "bb=" << endl; - octave_print_internal(cout,bb,0); - if(compz == 'V') - { - cout << endl << "ZZ=" << endl; - octave_print_internal(cout,ZZ,0); - } - cout << endl; - #endif + switch (ord_job[0]) + { + case 'S': + sort_test = &fin; + break; - // manually update alphar, alphai, betar - static int jj; - jj=0; - while(jj < nn) - { - #ifdef DEBUG_EIG - cout << "computing gen eig #" << jj << endl; - #endif - - static int zcnt; // number of zeros in this block - if(jj == (nn-1)) - zcnt = 1; - else if(aa(jj+1,jj) == 0) - zcnt = 1; - else zcnt = 2; - - if(zcnt == 1) // real zero - { - #ifdef DEBUG_EIG - cout << " single gen eig:" << endl; - cout << " alphar(" << jj << ") = " << aa(jj,jj) << endl; - cout << " betar( " << jj << ") = " << bb(jj,jj) << endl; - cout << " alphai(" << jj << ") = 0" << endl; - #endif + case 'B': + sort_test = &fout; + break; - alphar(jj) = aa(jj,jj); - alphai(jj) = 0; - betar(jj) = bb(jj,jj); - } - else // complex conjugate pair - { - #ifdef DEBUG_EIG - cout << "qz: calling dlag2:" << endl; - cout << "safmin=" << setiosflags(ios::scientific) << safmin << endl; - for(int idr = jj ; idr <= jj+1 ; idr++) - { - for(int idc = jj ; idc <= jj+1 ; idc++) - { - cout << "aa(" << idr << "," << idc << ")=" - << aa(idr,idc) << endl; - cout << "bb(" << idr << "," << idc << ")=" - << bb(idr,idc) << endl; - } - } - #endif - double scale1, scale2, wr1, wr2, wi; - F77_XFCN( dlag2, DLAG2, ( &aa(jj,jj), nn, &bb(jj,jj), nn, safmin, - scale1, scale2, wr1, wr2, wi)); + case '+': + sort_test = &fcrhp; + break; - #ifdef DEBUG_EIG - cout << "dlag2 returns: scale1=" << scale1 - << "\tscale2=" << scale2 << endl - << "\twr1=" << wr1 << "\twr2=" << wr2 - << "\twi=" << wi << endl; - #endif - // just to be safe, check if it's a real pair - if(wi == 0) - { - alphar(jj) = wr1; - alphai(jj) = 0; - betar(jj) = scale1; - alphar(jj+1) = wr2; - alphai(jj+1) = 0; - betar(jj+1) = scale2; - } - else - { - alphar(jj) = alphar(jj+1)=wr1; - alphai(jj) = -(alphai(jj+1) = wi); - betar(jj) = betar(jj+1) = scale1; - } - } - - jj += zcnt; // advance past this block + case '-': + sort_test = &folhp; + break; + + default: + // illegal order option (should never happen, since we + // checked the options at the top). + panic_impossible (); + break; } - #ifdef DEBUG_SORT - cout << "qz: back from dsubsp: aa=" << endl; - octave_print_internal(cout,aa,0); - cout << endl << "bb=" << endl; - octave_print_internal(cout,bb,0); - if(compz == 'V') - { - cout << endl << "ZZ=" << endl; - octave_print_internal(cout,ZZ,0); - } - cout << endl << "qz: ndim=" << ndim << endl << "fail=" << fail << endl; - cout << "alphar = " << endl; - octave_print_internal(cout,(Matrix) alphar,0); - cout << endl << "alphai = " << endl; - octave_print_internal(cout,(Matrix) alphai,0); - cout << endl << "beta = " << endl; - octave_print_internal(cout,(Matrix) betar,0); - cout << endl; - #endif + int ndim, fail, ind[nn]; + double inf_norm; + + F77_XFCN (xdlange, XDLANGE, + ("I", nn, nn, aa.fortran_vec (), nn, + work.fortran_vec (), inf_norm)); + + double eps = DBL_EPSILON*inf_norm*nn; + +#ifdef DEBUG_SORT + cout << "qz: calling dsubsp: aa=" << endl; + octave_print_internal (cout, aa, 0); + cout << endl << "bb=" << endl; + octave_print_internal (cout, bb, 0); + if (compz == 'V') + { + cout << endl << "ZZ=" << endl; + octave_print_internal (cout, ZZ, 0); + } + cout << endl; + cout << "alphar = " << endl; + octave_print_internal (cout, (Matrix) alphar, 0); + cout << endl << "alphai = " << endl; + octave_print_internal (cout, (Matrix) alphai, 0); + cout << endl << "beta = " << endl; + octave_print_internal (cout, (Matrix) betar, 0); + cout << endl; +#endif + + F77_XFCN (dsubsp, DSUBSP, + (nn, nn, aa.fortran_vec(), bb.fortran_vec(), + ZZ.fortran_vec(), sort_test, eps, ndim, fail, ind)); + +#ifdef DEBUG + cout << "qz: back from dsubsp: aa=" << endl; + octave_print_internal (cout, aa, 0); + cout << endl << "bb=" << endl; + octave_print_internal (cout, bb, 0); + if (compz == 'V') + { + cout << endl << "ZZ=" << endl; + octave_print_internal (cout, ZZ, 0); + } + cout << endl; +#endif + + // manually update alphar, alphai, betar + static int jj; + + jj=0; + while (jj < nn) + { +#ifdef DEBUG_EIG + cout << "computing gen eig #" << jj << endl; +#endif + + static int zcnt; // number of zeros in this block + + if (jj == (nn-1)) + zcnt = 1; + else if (aa(jj+1,jj) == 0) + zcnt = 1; + else zcnt = 2; + + if (zcnt == 1) // real zero + { +#ifdef DEBUG_EIG + cout << " single gen eig:" << endl; + cout << " alphar(" << jj << ") = " << aa(jj,jj) << endl; + cout << " betar( " << jj << ") = " << bb(jj,jj) << endl; + cout << " alphai(" << jj << ") = 0" << endl; +#endif + + alphar(jj) = aa(jj,jj); + alphai(jj) = 0; + betar(jj) = bb(jj,jj); + } + else + { + // complex conjugate pair +#ifdef DEBUG_EIG + cout << "qz: calling dlag2:" << endl; + cout << "safmin=" + << setiosflags (ios::scientific) << safmin << endl; + + for (int idr = jj; idr <= jj+1; idr++) + { + for (int idc = jj; idc <= jj+1; idc++) + { + cout << "aa(" << idr << "," << idc << ")=" + << aa(idr,idc) << endl; + cout << "bb(" << idr << "," << idc << ")=" + << bb(idr,idc) << endl; + } + } +#endif + + double scale1, scale2, wr1, wr2, wi; + F77_XFCN (dlag2, DLAG2, + (&aa(jj,jj), nn, &bb(jj,jj), nn, safmin, + scale1, scale2, wr1, wr2, wi)); + +#ifdef DEBUG_EIG + cout << "dlag2 returns: scale1=" << scale1 + << "\tscale2=" << scale2 << endl + << "\twr1=" << wr1 << "\twr2=" << wr2 + << "\twi=" << wi << endl; +#endif + + // just to be safe, check if it's a real pair + if (wi == 0) + { + alphar(jj) = wr1; + alphai(jj) = 0; + betar(jj) = scale1; + alphar(jj+1) = wr2; + alphai(jj+1) = 0; + betar(jj+1) = scale2; + } + else + { + alphar(jj) = alphar(jj+1)=wr1; + alphai(jj) = -(alphai(jj+1) = wi); + betar(jj) = betar(jj+1) = scale1; + } + } + + // advance past this block + jj += zcnt; + } + +#ifdef DEBUG_SORT + cout << "qz: back from dsubsp: aa=" << endl; + octave_print_internal (cout, aa, 0); + cout << endl << "bb=" << endl; + octave_print_internal (cout, bb, 0); + + if (compz == 'V') + { + cout << endl << "ZZ=" << endl; + octave_print_internal (cout, ZZ, 0); + } + cout << endl << "qz: ndim=" << ndim << endl + << "fail=" << fail << endl; + cout << "alphar = " << endl; + octave_print_internal (cout, (Matrix) alphar, 0); + cout << endl << "alphai = " << endl; + octave_print_internal (cout, (Matrix) alphai, 0); + cout << endl << "beta = " << endl; + octave_print_internal (cout, (Matrix) betar, 0); + cout << endl; +#endif + } } - } - + // compute generalized eigenvalues? ComplexColumnVector gev; - if(nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4)) - { - if(complex_case) - error("complex case not yet implemented"); - else + + if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4)) { - #ifdef DEBUG - cout << "qz: computing generalized eigenvalues" << endl; - #endif + if (complex_case) + { + error ("complex case not yet implemented"); + return retval; + } + else + { +#ifdef DEBUG + cout << "qz: computing generalized eigenvalues" << endl; +#endif - // return finite generalized eigenvalues - int ii, cnt = 0; - for( ii=0 ; ii < nn ; ii++) - if(betar(ii) != 0) - cnt++; - ComplexColumnVector tmp(cnt); - for( ii=0 ; ii < nn ; ii++) - if(betar(ii) != 0) - tmp(ii) = Complex(alphar(ii), alphai(ii))/betar(ii); - gev = tmp; + // return finite generalized eigenvalues + int cnt = 0; + + for (int ii = 0; ii < nn; ii++) + if (betar(ii) != 0) + cnt++; + + ComplexColumnVector tmp(cnt); + + for (int ii = 0; ii < nn; ii++) + if (betar(ii) != 0) + tmp(ii) = Complex(alphar(ii), alphai(ii))/betar(ii); + gev = tmp; + } } - } // right, left eigenvector matrices - if(nargout >= 5) - { - char side = (nargout == 5 ? 'R' : 'B'), // which side to compute? - howmny = 'B'; // compute all of them and backtransform - int *select = NULL; // dummy pointer; select is not used. - int m; - - if(complex_case) - error("complex type not yet implemented"); - else + if (nargout >= 5) { - #ifdef DEBUG - cout << "qz: computing generalized eigenvectors" << endl; - #endif + char side = (nargout == 5 ? 'R' : 'B'); // which side to compute? + char howmny = 'B'; // compute all of them and backtransform + int *select = NULL; // dummy pointer; select is not used. + int m; + + if (complex_case) + { + error ("complex type not yet implemented"); + return retval; + } + else + { +#ifdef DEBUG + cout << "qz: computing generalized eigenvectors" << endl; +#endif + + VL = QQ; + VR = ZZ; + + F77_XFCN (dtgevc, DTGEVC, + (&side, &howmny, select, nn, aa.fortran_vec(), + nn, bb.fortran_vec(), nn, VL.fortran_vec(), nn, + VR.fortran_vec(), nn, nn, m, work.fortran_vec(), + info, 1L, 1L)); + + if (f77_exception_encountered) + { + error ("unrecoverable error in qz (dtgevc)"); + return retval; + } - VL = QQ; - VR = ZZ; - - F77_XFCN( dtgevc, DTGEVC, ( &side, &howmny, select, nn, aa.fortran_vec(), - nn, bb.fortran_vec(), nn, VL.fortran_vec(), nn, VR.fortran_vec(), - nn, nn, m, work.fortran_vec(), info, 1L, 1L )); - if(f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in qz(dtgevc)"); + // now construct the complex form of VV, WW + int jj = 0; + + while (jj < nn) + { + // see if real or complex eigenvalue + int cinc = 2; // column increment; assume complex eigenvalue + + if (jj == (nn-1)) + cinc = 1; // single column + else if (aa(jj+1,jj) == 0) + cinc = 1; + + // now copy the eigenvector (s) to CVR, CVL + if (cinc == 1) + { + for (int ii = 0; ii < nn; ii++) + CVR(ii,jj) = VR(ii,jj); + + if (side == 'B') + for (int ii = 0; ii < nn; ii++) + CVL(ii,jj) = VL(ii,jj); + } + else + { + // double column; complex vector + + for (int ii = 0; ii < nn; ii++) + { + CVR(ii,jj) = Complex (VR(ii,jj), VR(ii,jj+1)); + CVR(ii,jj+1) = Complex (VR(ii,jj), -VR(ii,jj+1)); + } - // now construct the complex form of VV, WW - int jj = 0; - while(jj < nn) - { - // see if real or complex eigenvalue - int cinc = 2; // column increment; assume complex eigenvalue - if(jj == (nn-1)) - cinc = 1; // single column - else if(aa(jj+1,jj) == 0) - cinc = 1; - - // now copy the eigenvector(s) to CVR, CVL - if(cinc == 1) - { - int ii; - for(ii = 0; ii < nn ; ii++) - CVR(ii,jj) = VR(ii,jj); - if(side == 'B') - for(ii = 0; ii < nn ; ii++) - CVL(ii,jj) = VL(ii,jj); - } - else // double column; complex vector - { - int ii; - for(ii = 0; ii < nn ; ii++) - { - CVR(ii,jj) = Complex(VR(ii,jj),VR(ii,jj+1)); - CVR(ii,jj+1) = Complex(VR(ii,jj),-VR(ii,jj+1)); - } - if(side == 'B') - for(ii = 0; ii < nn ; ii++) - { - CVL(ii,jj) = Complex(VL(ii,jj),VL(ii,jj+1)); - CVL(ii,jj+1) = Complex(VL(ii,jj),-VL(ii,jj+1)); - } - } - jj += cinc; // advance to next eigenvectors (if any) - } + if (side == 'B') + for (int ii = 0; ii < nn; ii++) + { + CVL(ii,jj) = Complex (VL(ii,jj), VL(ii,jj+1)); + CVL(ii,jj+1) = Complex (VL(ii,jj), -VL(ii,jj+1)); + } + } + + // advance to next eigenvectors (if any) + jj += cinc; + } + } } + + switch (nargout) + { + case 7: + retval(6) = gev; + + case 6: // return eigenvectors + retval(5) = CVL; + + case 5: // return eigenvectors + retval(4) = CVR; + + case 4: + if (nargin == 3) + { +#ifdef DEBUG + cout << "qz: sort: retval(3) = gev = " << endl; + octave_print_internal (cout, gev); + cout << endl; +#endif + retval(3) = gev; + } + else + retval(3) = ZZ; + + case 3: + if (nargin == 3) + retval(2) = ZZ; + else + retval(2) = QQ; + + case 2: +#ifdef DEBUG + cout << "qz: retval (1) = bb = " << endl; + octave_print_internal (cout, bb, 0); + cout << endl << "qz: retval(0) = aa = " <<endl; + octave_print_internal (cout, aa, 0); + cout << endl; +#endif + retval(1) = bb; + retval(0) = aa; + break; + + case 1: + case 0: +#ifdef DEBUG + cout << "qz: retval(0) = gev = " << gev << endl; +#endif + retval(0) = gev; + break; + + default: + error ("qz: too many return arguments."); + break; } - switch(nargout) - { - case 7: - retval(6) = gev; - case 6: // return eigenvectors - retval(5) = CVL; - case 5: // return eigenvectors - retval(4) = CVR; - case 4: - if(nargin == 3) - { - #ifdef DEBUG - cout << "qz: sort: retval(3) = gev = " << endl; - octave_print_internal(cout,gev); - cout << endl; - #endif - retval(3) = gev; - } - else retval(3) = ZZ; - case 3: - if(nargin == 3) - retval(2) = ZZ; - else - retval(2) = QQ; - case 2: - #ifdef DEBUG - cout << "qz: retval(1) = bb = " << endl; - octave_print_internal(cout,bb,0); - cout << endl << "qz: retval(0) = aa = " <<endl; - octave_print_internal(cout,aa,0); - cout << endl; - #endif - retval(1) = bb; - retval(0) = aa; - break; - case 1: - case 0: - #ifdef DEBUG - cout << "qz: retval(0) = gev = " << gev << endl; - #endif - retval(0) = gev; - break; - default: - error("qz: too many return arguments. Sorry. "); - } +#ifdef DEBUG + cout << "qz: exiting (at long last)" << endl; +#endif - #ifdef DEBUG - cout << "qz: exiting (at long last)" << endl; - #endif - return retval; }
--- a/src/Makefile.in +++ b/src/Makefile.in @@ -43,7 +43,7 @@ det.cc eig.cc expm.cc fft.cc fft2.cc filter.cc find.cc \ fsolve.cc gammainc.cc getgrent.cc getpwent.cc getrusage.cc \ givens.cc hess.cc ifft.cc ifft2.cc inv.cc log.cc lpsolve.cc \ - lsode.cc lu.cc minmax.cc pinv.cc qr.cc quad.cc qzval.cc rand.cc \ + lsode.cc lu.cc minmax.cc pinv.cc qr.cc quad.cc qz.cc rand.cc \ schur.cc sort.cc svd.cc syl.cc time.cc DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))
--- a/src/defaults.cc +++ b/src/defaults.cc @@ -46,6 +46,7 @@ #include "file-ops.h" #include "gripes.h" #include "help.h" +#include "oct-obj.h" #include "ov.h" #include "toplev.h" #include "variables.h" @@ -427,6 +428,16 @@ "Octave version"); } +DEFUN (rehash, , , + "rehash (): reinitialize LOADPATH directory cache") +{ + octave_value_list retval; + + Vload_path_dir_path.rehash (); + + return retval; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/src/help.cc +++ b/src/help.cc @@ -640,7 +640,8 @@ } else { - error ("help: unable to find info!"); + error ("help: unable to find info"); + error ("help: you need info 2.18 or later (texinfo 3.12)"); break; } }
--- a/src/load-save.cc +++ b/src/load-save.cc @@ -27,6 +27,7 @@ #include <cfloat> #include <cstring> #include <cctype> +#include <ctime> #include <string> @@ -41,6 +42,7 @@ #include "glob-match.h" #include "lo-mappers.h" #include "mach-info.h" +#include "oct-env.h" #include "str-vec.h" #include "defun.h" @@ -55,6 +57,7 @@ #include "unwind-prot.h" #include "utils.h" #include "variables.h" +#include "version.h" // The default output format. May be one of "binary", "text", or // "mat-binary". @@ -2366,28 +2369,53 @@ } static void -write_binary_header (ostream& os, load_save_format format) +write_header (ostream& os, load_save_format format) { - if (format == LS_BINARY) + switch (format) { - os << (oct_mach_info::words_big_endian () - ? "Octave-1-B" : "Octave-1-L"); - - oct_mach_info::float_format flt_fmt = - oct_mach_info::native_float_format (); - - char tmp = (char) float_format_to_mopt_digit (flt_fmt); - - os.write (&tmp, 1); + case LS_BINARY: + { + os << (oct_mach_info::words_big_endian () + ? "Octave-1-B" : "Octave-1-L"); + + oct_mach_info::float_format flt_fmt = + oct_mach_info::native_float_format (); + + char tmp = (char) float_format_to_mopt_digit (flt_fmt); + + os.write (&tmp, 1); + } + break; + + case LS_ASCII: + { + time_t now = time (0); + + string time_string = asctime (gmtime (&now)); + time_string = time_string.substr (0, time_string.length () - 1); + + os << "# Created by Octave " OCTAVE_VERSION ", " + << time_string + << " <" + << octave_env::get_user_name () + << "@" + << octave_env::get_host_name () + << ">" << "\n"; + } + break; + + default: + break; } } static void save_vars (const string_vector& argv, int argv_idx, int argc, ostream& os, bool save_builtins, load_save_format fmt, - bool save_as_floats) + bool save_as_floats, bool write_header_info) { - write_binary_header (os, fmt); + if (write_header_info) + write_header (os, fmt); if (argv_idx == argc) { @@ -2424,7 +2452,7 @@ if (file) { - save_vars (string_vector (), 0, 0, file, false, format, false); + save_vars (string_vector (), 0, 0, file, false, format, false, true); message (0, "save to `%s' complete", fname); } else @@ -2432,7 +2460,7 @@ } DEFUN_TEXT (save, args, , - "save [-ascii] [-binary] [-float-binary] [-mat-binary] \n\ + "save [-append] [-ascii] [-binary] [-float-binary] [-mat-binary] \n\ [-save-builtins] file [pattern ...]\n\ \n\ save variables in a file") @@ -2455,9 +2483,15 @@ load_save_format format = get_default_save_format (); + bool append = false; + int i; for (i = 1; i < argc; i++) { + if (argv[i] == "-append") + { + append = true; + } if (argv[i] == "-ascii" || argv[i] == "-a") { format = LS_ASCII; @@ -2503,7 +2537,7 @@ // in a octave_value (string)? save_vars (argv, i, argc, octave_stdout, save_builtins, format, - save_as_floats); + save_as_floats, true); } // Guard against things like `save a*', which are probably mistakes... @@ -2519,16 +2553,21 @@ i++; - unsigned mode = ios::out|ios::trunc; + unsigned mode = ios::out; if (format == LS_BINARY || format == LS_MAT_BINARY) mode |= ios::bin; + mode |= append ? ios::ate : ios::trunc; + ofstream file (fname.c_str (), mode); if (file) { + bool write_header_info + = ((file.rdbuf ())->seekoff (0, ios::cur) == 0); + save_vars (argv, i, argc, file, save_builtins, format, - save_as_floats); + save_as_floats, write_header_info); } else {