Mercurial > hg > octave-nkf
diff liboctave/dMatrix.cc @ 3130:02766207b74c
[project @ 1998-01-25 08:27:23 by jwe]
author | jwe |
---|---|
date | Sun, 25 Jan 1998 08:27:25 +0000 |
parents | a6a00badcc12 |
children | 0d640dc625c7 |
line wrap: on
line diff
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -95,9 +95,9 @@ const double*, const int&, double&, int&, long, long); - double F77_FCN (dlange, DLANGE) (const char*, const int&, - const int&, const double*, - const int&, double*); + 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*); @@ -1334,18 +1334,22 @@ int nc = columns (); + // Preconditioning step 1: trace normalization to reduce dynamic + // range of poles, but avoid making stable eigenvalues unstable. + // trace shift value - double trshift = 0; - - // Preconditioning step 1: trace normalization. + double trshift = 0.0; for (int i = 0; i < nc; i++) trshift += m.elem (i, i); trshift /= nc; - for (int i = 0; i < nc; i++) - m.elem (i, i) -= trshift; + if (trshift > 0.0) + { + for (int i = 0; i < nc; i++) + m.elem (i, i) -= trshift; + } // Preconditioning step 2: balancing. @@ -1356,9 +1360,10 @@ // Preconditioning step 3: scaling. ColumnVector work(nc); - double inf_norm - = F77_FCN (dlange, DLANGE) ("I", nc, nc, m.fortran_vec (),nc, - work.fortran_vec ()); + double inf_norm; + + F77_FCN (xdlange, XDLANGE) ("I", nc, nc, m.fortran_vec (), nc, + work.fortran_vec (), inf_norm); int sqpow = (int) (inf_norm > 0.0 ? (1.0 + log (inf_norm) / log (2.0)) @@ -1396,7 +1401,7 @@ // Zero power. dpp = -dpp; - for(int j = 0; j < nc; j++) + for (int j = 0; j < nc; j++) { npp.elem (j, j) += 1.0; dpp.elem (j, j) += 1.0; @@ -1424,7 +1429,10 @@ // Reverse preconditioning step 1: fix trace normalization. - return retval * exp (trshift); + if (trshift > 0.0) + retval = exp (trshift) * retval; + + return retval; } Matrix&