Mercurial > hg > octave-terminal
changeset 22:2cd2476fb32d
[project @ 1993-08-10 20:28:05 by jwe]
Add init functions for balancing classes.
author | jwe |
---|---|
date | Tue, 10 Aug 1993 20:28:05 +0000 |
parents | 1b59f5c6c370 |
children | 36bbfe8c8d9f |
files | liboctave/Matrix-ext.cc |
diffstat | 1 files changed, 202 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/Matrix-ext.cc +++ b/liboctave/Matrix-ext.cc @@ -29,6 +29,206 @@ #include "mx-inlines.cc" /* + * AEPBALANCE operations + */ + +int +AEPBALANCE::init (const Matrix& a, const char *balance_job) +{ + if (a.nr != a.nc) + FAIL; + + int n = a.nc; + +// Parameters for balance call. + + int info; + int ilo; + int ihi; + double *scale = new double [n]; + +// Copy matrix into local structure. + + balanced_mat = a; + + F77_FCN (dgebal) (balance_job, &n, balanced_mat.fortran_vec (), + &n, &ilo, &ihi, scale, &info, 1L, 1L); + +// Initialize balancing matrix to identity. + + balancing_mat = Matrix (n, n, 0.0); + for (int i = 0; i < n; i++) + balancing_mat.elem (i ,i) = 1.0; + + F77_FCN (dgebak) (balance_job, "R", &n, &ilo, &ihi, scale, &n, + balancing_mat.fortran_vec (), &n, &info, 1L, 1L); + + delete [] scale; + + return info; +} + +int +ComplexAEPBALANCE::init (const ComplexMatrix& a, const char *balance_job) +{ + + int n = a.nc; + +// Parameters for balance call. + + int info; + int ilo; + int ihi; + double *scale = new double [n]; + +// Copy matrix into local structure. + + balanced_mat = a; + + F77_FCN (zgebal) (balance_job, &n, balanced_mat.fortran_vec (), + &n, &ilo, &ihi, scale, &info, 1L, 1L); + +// Initialize balancing matrix to identity. + + balancing_mat = Matrix (n, n, 0.0); + for (int i = 0; i < n; i++) + balancing_mat (i, i) = 1.0; + + F77_FCN (zgebak) (balance_job, "R", &n, &ilo, &ihi, scale, &n, + balancing_mat.fortran_vec(), &n, &info, 1L, 1L); + + delete [] scale; + + return info; +} + +/* + * GEPBALANCE operations + */ + +int +GEPBALANCE::init (const Matrix& a, const Matrix& b, const char *balance_job) +{ + if (a.nr != a.nc || a.nr != b.nr || b.nr != b.nc) + FAIL; + + int n = a.nc; + +// Parameters for balance call. + + int info; + int ilo; + int ihi; + double *cscale = new double [n]; + double *cperm = new double [n]; + Matrix wk (n, 6, 0.0); + +// Back out the permutations: +// +// cscale contains the exponents of the column scaling factors in its +// ilo through ihi locations and the reducing column permutations in +// its first ilo-1 and its ihi+1 through n locations. +// +// cperm contains the column permutations applied in grading the a and b +// submatrices in its ilo through ihi locations. +// +// wk contains the exponents of the row scaling factors in its ilo +// through ihi locations, the reducing row permutations in its first +// ilo-1 and its ihi+1 through n locations, and the row permutations +// applied in grading the a and b submatrices in its n+ilo through +// n+ihi locations. + +// Copy matrices into local structure. + + 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; + +// Check for permutation option. + + if (*balance_job == 'P' || *balance_job == 'B') + { + F77_FCN(reduce)(&n, &n, balanced_a_mat.fortran_vec (), + &n, balanced_b_mat.fortran_vec (), &ilo, &ihi, + cscale, wk.fortran_vec ()); + } + else + { + +// Set up for scaling later. + + ilo = 1; + ihi = n; + } + +// Check for scaling option. + + if ((*balance_job == 'S' || *balance_job == 'B') && ilo != ihi) + { + F77_FCN(scaleg)(&n, &n, balanced_a_mat.fortran_vec (), + &n, balanced_b_mat.fortran_vec (), &ilo, &ihi, + cscale, cperm, wk.fortran_vec ()); + } + else + { + +// Set scaling data to 0's. + + for (int tmp = ilo-1; tmp < ihi; tmp++) + { + cscale[tmp] = 0.0; + wk.elem(tmp,0) = 0.0; + } + } + +// 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) (balance_job, "R", &n, &ilo, &ihi, cscale, &n, + right_balancing_mat.fortran_vec (), &n, &info, 1L, + 1L); + +// Row permutations/scaling. + + F77_FCN (dgebak) (balance_job, "L", &n, &ilo, &ihi, &wk.elem (0, 0), &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 ((*balance_job == 'P' || *balance_job == 'B') && ilo != ihi) + { + F77_FCN (gradeq) (&n, &n, balanced_a_mat.fortran_vec (), + &n, balanced_b_mat.fortran_vec (), &ilo, &ihi, + cperm, &wk.elem (0, 1)); + } +#endif + +// Transpose for aa = cc*a*dd convention... + left_balancing_mat = left_balancing_mat.transpose (); + + delete [] cscale; + delete [] cperm; + + return info; +} + +/* * HESS stuff */ @@ -38,7 +238,7 @@ if (a.nr != a.nc) FAIL; - char jobbal = 'S'; + char jobbal = 'N'; char side = 'R'; int n = a.nc; @@ -96,7 +296,7 @@ if (a.nr != a.nc) FAIL; - char job = 'S'; + char job = 'N'; char side = 'R'; int n = a.nc;