diff liboctave/CMatrix.cc @ 3130:02766207b74c

[project @ 1998-01-25 08:27:23 by jwe]
author jwe
date Sun, 25 Jan 1998 08:27:25 +0000
parents 528f4270e904
children fccab8e7d35f
line wrap: on
line diff
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -100,9 +100,9 @@
 				const Complex*, const int&, double&,
 				int&, long, long);
 
-  double F77_FCN (zlange, ZLANGE) (const char*, const int&,
-				   const int&, const Complex*,
-				   const int&, double*); 
+  int F77_FCN (xzlange, XZLANGE) (const char*, const int&,
+				  const int&, const Complex*,
+				  const int&, double*, double&); 
 }
 
 static const Complex Complex_NaN_result (octave_NaN, octave_NaN);
@@ -1586,16 +1586,20 @@
 
   int nc = columns ();
 
+  // Preconditioning step 1: trace normalization to reduce dynamic
+  // range of poles, but avoid making stable eigenvalues unstable.
+
   // trace shift value
   Complex trshift = 0.0;
 
-  // Preconditioning step 1: trace normalization.
-
   for (int i = 0; i < nc; i++)
     trshift += m.elem (i, i);
 
   trshift /= nc;
 
+  if (trshift.real () < 0.0)
+    trshift = trshift.imag ();
+
   for (int i = 0; i < nc; i++)
     m.elem (i, i) -= trshift;
 
@@ -1608,9 +1612,10 @@
   // Preconditioning step 3: scaling.
 
   ColumnVector work (nc);
-  double inf_norm
-    = F77_FCN (zlange, ZLANGE) ("I", nc, nc, m.fortran_vec (), nc,
-				work.fortran_vec ());
+  double inf_norm;
+
+  F77_FCN (xzlange, XZLANGE) ("I", nc, nc, m.fortran_vec (), nc,
+			      work.fortran_vec (), inf_norm);
 
   int sqpow = (inf_norm > 0.0
 	       ? static_cast<int> (1.0 + log (inf_norm) / log (2.0)) : 0);
@@ -1677,7 +1682,7 @@
 
   // Reverse preconditioning step 1: fix trace normalization.
 
-  return retval * exp (trshift);
+  return exp (trshift) * retval;
 }
 
 // column vector by row vector -> matrix operations