Mercurial > octave
diff src/DLD-FUNCTIONS/balance.cc @ 3185:9580887dd160
[project @ 1998-09-26 02:45:55 by jwe]
author | jwe |
---|---|
date | Sat, 26 Sep 1998 02:45:59 +0000 |
parents | 3f6a813eb09e |
children | eba59b8c64dc |
line wrap: on
line diff
--- a/src/DLD-FUNCTIONS/balance.cc Fri Sep 25 04:35:09 1998 +0000 +++ b/src/DLD-FUNCTIONS/balance.cc Sat Sep 26 02:45:59 1998 +0000 @@ -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; }