Mercurial > hg > octave-lojdl
changeset 1933:333923894848
[project @ 1996-02-12 03:07:39 by jwe]
author | jwe |
---|---|
date | Mon, 12 Feb 1996 03:07:39 +0000 |
parents | 682f31b20894 |
children | 0e591d443ff0 |
files | liboctave/CmplxAEPBAL.cc liboctave/CmplxAEPBAL.h liboctave/dbleAEPBAL.cc liboctave/dbleAEPBAL.h liboctave/dbleGEPBAL.cc liboctave/dbleGEPBAL.h |
diffstat | 6 files changed, 156 insertions(+), 112 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CmplxAEPBAL.cc +++ b/liboctave/CmplxAEPBAL.cc @@ -50,44 +50,47 @@ int ComplexAEPBALANCE::init (const ComplexMatrix& a, const string& balance_job) { - int a_nc = a.cols (); + int n = a.cols (); - if (a.rows () != a_nc) + if (a.rows () != n) { (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix"); return -1; } - int n = a_nc; - - // Parameters for balance call. - int info; int ilo; int ihi; - double *scale = new double [n]; - // Copy matrix into local structure. + Array<double> scale (n); + double *pscale = scale.fortran_vec (); balanced_mat = a; - - char bal_job = balance_job[0]; + Complex *p_balanced_mat = balanced_mat.fortran_vec (); - F77_FCN (zgebal, ZGEBAL) (&bal_job, n, - balanced_mat.fortran_vec (), n, ilo, ihi, - scale, info, 1L, 1L); + char job = balance_job[0]; - // Initialize balancing matrix to identity. + F77_XFCN (zgebal, ZGEBAL, (&job, n, p_balanced_mat, n, ilo, ihi, + pscale, info, 1L, 1L)); - balancing_mat = Matrix (n, n, 0.0); - for (int i = 0; i < n; i++) - balancing_mat (i, i) = 1.0; + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in zgebal"); + else + { + balancing_mat = Matrix (n, n, 0.0); + for (int i = 0; i < n; i++) + balancing_mat.elem (i, i) = 1.0; - F77_FCN (zgebak, ZGEBAK) (&bal_job, "R", n, ilo, ihi, scale, n, - balancing_mat.fortran_vec (), n, info, 1L, - 1L); + Complex *p_balancing_mat = balancing_mat.fortran_vec (); + + char side = 'R'; - delete [] scale; + F77_XFCN (zgebak, ZGEBAK, (&job, &side, n, ilo, ihi, pscale, n, + p_balancing_mat, n, info, 1L, 1L)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in zgebak"); + } return info; }
--- a/liboctave/CmplxAEPBAL.h +++ b/liboctave/CmplxAEPBAL.h @@ -59,6 +59,8 @@ return *this; } + ~ComplexAEPBALANCE (void) { } + ComplexMatrix balanced_matrix (void) const { return balanced_mat; } ComplexMatrix balancing_matrix (void) const { return balancing_mat; }
--- a/liboctave/dbleAEPBAL.cc +++ b/liboctave/dbleAEPBAL.cc @@ -49,44 +49,47 @@ int AEPBALANCE::init (const Matrix& a, const string& balance_job) { - int a_nc = a.cols (); + int n = a.cols (); - if (a.rows () != a_nc) + if (a.rows () != n) { (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix"); return -1; } - int n = a_nc; - - // Parameters for balance call. - int info; int ilo; int ihi; - double *scale = new double [n]; - // Copy matrix into local structure. + Array<double> scale (n); + double *pscale = scale.fortran_vec (); balanced_mat = a; - - char bal_job = balance_job[0]; + double *p_balanced_mat = balanced_mat.fortran_vec (); - F77_FCN (dgebal, DGEBAL) (&bal_job, n, - balanced_mat.fortran_vec (), n, ilo, ihi, - scale, info, 1L, 1L); + char job = balance_job[0]; - // Initialize balancing matrix to identity. + F77_XFCN (dgebal, DGEBAL, (&job, n, p_balanced_mat, n, ilo, ihi, + pscale, info, 1L, 1L)); - balancing_mat = Matrix (n, n, 0.0); - for (int i = 0; i < n; i++) - balancing_mat.elem (i ,i) = 1.0; + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in dgebal"); + else + { + balancing_mat = Matrix (n, n, 0.0); + for (int i = 0; i < n; i++) + balancing_mat.elem (i ,i) = 1.0; - F77_FCN (dgebak, DGEBAK) (&bal_job, "R", n, ilo, ihi, scale, n, - balancing_mat.fortran_vec (), n, info, 1L, - 1L); + double *p_balancing_mat = balancing_mat.fortran_vec (); + + char side = 'R'; - delete [] scale; + F77_XFCN (dgebak, DGEBAK, (&job, &side, n, ilo, ihi, pscale, n, + p_balancing_mat, n, info, 1L, 1L)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in dgebak"); + } return info; }
--- a/liboctave/dbleAEPBAL.h +++ b/liboctave/dbleAEPBAL.h @@ -59,6 +59,8 @@ return *this; } + ~AEPBALANCE (void) { } + Matrix balanced_matrix (void) const { return balanced_mat; } Matrix balancing_matrix (void) const { return balancing_mat; }
--- a/liboctave/dbleGEPBAL.cc +++ b/liboctave/dbleGEPBAL.cc @@ -73,14 +73,15 @@ int n = a_nc; - // Parameters for balance call. - int info; int ilo; int ihi; - double *cscale = new double [n]; - double *cperm = new double [n]; + + Array<double> cscale (n); + double *pcscale = cscale.fortran_vec (); + Matrix wk (n, 6, 0.0); + double *pwk = wk.fortran_vec (); // Back out the permutations: // @@ -102,23 +103,18 @@ balanced_a_mat = a; balanced_b_mat = b; - // Initialize balancing matrices to identity. - - left_balancing_mat = Matrix (n, n, 0.0); - for (int i = 0; i < n; i++) - left_balancing_mat (i, i) = 1.0; - - right_balancing_mat = left_balancing_mat; + double *p_balanced_a_mat = balanced_a_mat.fortran_vec (); + double *p_balanced_b_mat = balanced_b_mat.fortran_vec (); // Check for permutation option. - char bal_job = balance_job[0]; + char job = balance_job[0]; - if (bal_job == 'P' || bal_job == 'B') + if (job == 'P' || job == 'B') { - F77_FCN (reduce, REDUCE) (n, n, balanced_a_mat.fortran_vec (), - n, balanced_b_mat.fortran_vec (), ilo, - ihi, cscale, wk.fortran_vec ()); + F77_XFCN (reduce, REDUCE, (n, n, p_balanced_a_mat, n, + p_balanced_b_mat, ilo, ihi, + pcscale, pwk)); } else { @@ -128,66 +124,102 @@ ihi = n; } - // Check for scaling option. - - if ((bal_job == 'S' || bal_job == 'B') && ilo != ihi) - { - F77_FCN (scaleg, SCALEG) (n, n, balanced_a_mat.fortran_vec (), - n, balanced_b_mat.fortran_vec (), ilo, - ihi, cscale, cperm, wk.fortran_vec ()); - } + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in reduce"); else { - // Set scaling data to 0's. + Array<double> cperm (n); + double *pcperm = cperm.fortran_vec (); + + // Check for scaling option. + + if ((job == 'S' || job == 'B') && ilo != ihi) + { + F77_XFCN (scaleg, SCALEG, (n, n, p_balanced_a_mat, n, + p_balanced_b_mat, ilo, ihi, + pcscale, pcperm, pwk)); + } + else + { + // Set scaling data to 0's. + + for (int i = ilo-1; i < ihi; i++) + { + cscale.elem (i) = 0.0; + wk.elem (i, 0) = 0.0; + } + } + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in scaleg"); + else + { + // Scaleg returns exponents, not values, so... + + for (int i = ilo-1; i < ihi; i++) + { + cscale.elem (i) = pow (2.0, cscale.elem (i)); + wk.elem (i, 0) = pow (2.0, -wk.elem (i, 0)); + } + + // Initialize balancing matrices to identity. + + left_balancing_mat = Matrix (n, n, 0.0); + for (int i = 0; i < n; i++) + left_balancing_mat (i, i) = 1.0; + + right_balancing_mat = left_balancing_mat; + + double *p_right_balancing_mat = right_balancing_mat.fortran_vec (); + double *p_left_balancing_mat = left_balancing_mat.fortran_vec (); - for (int tmp = ilo-1; tmp < ihi; tmp++) - { - cscale[tmp] = 0.0; - wk.elem (tmp, 0) = 0.0; + // Column permutations/scaling. + + char side = 'R'; + + F77_XFCN (dgebak, DGEBAK, (&job, &side, n, ilo, ihi, pcscale, + n, p_right_balancing_mat, n, info, + 1L, 1L)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in dgebak"); + else + { + // Row permutations/scaling. + + side = 'L'; + + F77_XFCN (dgebak, DGEBAK, (&job, &side, n, ilo, ihi, pwk, + n, p_left_balancing_mat, n, + info, 1L, 1L)); + +#if 0 + // XXX FIXME XXX --- these four lines need to be added and + // debugged. GEPBALANCE::init will work without them, though, so + // here they are. + + if ((job == 'P' || job == 'B') && ilo != ihi) + { + F77_XFCN (gradeq, GRADEQ, (n, n, p_balanced_a_mat, n, + p_balanced_b_mat, ilo, ihi, + pcperm, pwk)); + } +#endif + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in dgebak"); + else + { + // Transpose for aa = cc*a*dd convention... + + left_balancing_mat = left_balancing_mat.transpose (); + } + } } } - // Scaleg returns exponents, not values, so... - - for (int tmp = ilo-1; tmp < ihi; tmp++) - { - cscale[tmp] = pow (2.0, cscale[tmp]); - wk.elem (tmp, 0) = pow (2.0, -wk.elem (tmp, 0)); - } - - // Column permutations/scaling. - - F77_FCN (dgebak, DGEBAK) (&bal_job, "R", n, ilo, ihi, cscale, n, - right_balancing_mat.fortran_vec (), n, - info, 1L, 1L); - - // Row permutations/scaling. - - F77_FCN (dgebak, DGEBAK) (&bal_job, "L", n, ilo, ihi, - wk.fortran_vec (), n, - left_balancing_mat.fortran_vec (), n, - info, 1L, 1L); - - // XXX FIXME XXX --- these four lines need to be added and - // debugged. GEPBALANCE::init will work without them, though, so - // here they are. - -#if 0 - if ((bal_job == 'P' || bal_job == 'B') && ilo != ihi) - { - F77_FCN (gradeq, GRADEQ) (n, n, balanced_a_mat.fortran_vec (), - n, balanced_b_mat.fortran_vec (), ilo, - ihi, cperm, wk.fortran_vec ()); - } -#endif - - // Transpose for aa = cc*a*dd convention... - - left_balancing_mat = left_balancing_mat.transpose (); - - delete [] cscale; - delete [] cperm; - return info; }