# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1223800274 -7200
# Node ID f53f57d2ee51e0de22bb8c214b8b7848cc7e789a
# Parent  757ca614e9181578be235ad677cc1c74a2dc03cc
improve inverse preconditioning according to Marco Caliari

diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -2964,21 +2964,30 @@
   // inverse scaling (diagonal transformation)
   for (octave_idx_type i = 0; i < nc; i++)
     for (octave_idx_type j = 0; j < nc; j++)
-       retval(i,j) *= dscale(i) / dscale(j);
+      retval(i,j) *= dscale(i) / dscale(j);
 
   OCTAVE_QUIT;
 
   // construct balancing permutation vector
   Array<octave_idx_type> iperm (nc);
   for (octave_idx_type i = 0; i < nc; i++)
-    iperm(i) = i;  // initialize to identity permutation
+    iperm(i) = i;  // identity permutation
+
+  // trailing permutations must be done in reverse order
+  for (octave_idx_type i = nc - 1; i >= ihi; i--)
+    {
+      octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
+      octave_idx_type tmp = iperm(i);
+      iperm(i) = iperm(swapidx);
+      iperm(swapidx) = tmp;
+    }
 
   // leading permutations in forward order
   for (octave_idx_type i = 0; i < (ilo-1); i++)
     {
       octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
       octave_idx_type tmp = iperm(i);
-      iperm(i) = iperm(swapidx);
+      iperm(i) = iperm (swapidx);
       iperm(swapidx) = tmp;
     }
 
@@ -2994,32 +3003,7 @@
     for (octave_idx_type j = 0; j < nc; j++)
       retval(i,j) = tmpMat(invpvec(i),invpvec(j));
 
-  OCTAVE_QUIT;
-
-  for (octave_idx_type i = 0; i < nc; i++)
-    iperm(i) = i;  // initialize to identity permutation
-
-  // trailing permutations must be done in reverse order
-  for (octave_idx_type i = nc - 1; i >= ihi; i--)
-    {
-      octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
-      octave_idx_type tmp = iperm(i);
-      iperm(i) = iperm(swapidx);
-      iperm(swapidx) = tmp;
-    }
-
-  // construct inverse balancing permutation vector
-  for (octave_idx_type i = 0; i < nc; i++)
-    invpvec(iperm(i)) = i;     // Thanks to R. A. Lippert for this method
-
-  OCTAVE_QUIT;
-
-  tmpMat = retval;
-  for (octave_idx_type i = 0; i < nc; i++)
-    for (octave_idx_type j = 0; j < nc; j++)
-      retval(i,j) = tmpMat(invpvec(i),invpvec(j));
-
-  // Reverse preconditioning step 1: fix trace normalization.
+  // Reverse preconditioning step 1: fix trace normalization. 
 
   return exp (trshift) * retval;
 }
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,11 @@
+2008-10-12  Jaroslav Hajek <highegg@gmail.com>
+
+	* CSparse.cc (ComplexMatrix::expm): Improve inverse preconditioning
+	according to Marco Caliari.
+	* dSparse.cc (Matrix::expm): Likewise.
+	* fCSparse.cc (FloatComplexMatrix::expm): Likewise.
+	* fSparse.cc (FloatMatrix::expm): Likewise.
+
 2008-10-10  Jaroslav Hajek  <highegg@gmail.com>
 
 	* sparse-util.h (SparseCholPrint): Change char * argument to const
diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -2589,7 +2589,7 @@
   // apply inverse scaling to computed exponential
   for (octave_idx_type i = 0; i < nc; i++)
     for (octave_idx_type j = 0; j < nc; j++)
-       retval(i,j) *= dscale(i) / dscale(j);
+      retval(i,j) *= dscale(i) / dscale(j);
 
   OCTAVE_QUIT;
 
@@ -2598,6 +2598,15 @@
   for (octave_idx_type i = 0; i < nc; i++)
     iperm(i) = i;  // identity permutation
 
+  // trailing permutations must be done in reverse order
+  for (octave_idx_type i = nc - 1; i >= ihi; i--)
+    {
+      octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
+      octave_idx_type tmp = iperm(i);
+      iperm(i) = iperm(swapidx);
+      iperm(swapidx) = tmp;
+    }
+
   // leading permutations in forward order
   for (octave_idx_type i = 0; i < (ilo-1); i++)
     {
@@ -2613,39 +2622,13 @@
     invpvec(iperm(i)) = i;     // Thanks to R. A. Lippert for this method
 
   OCTAVE_QUIT;
- 
+
   Matrix tmpMat = retval;
   for (octave_idx_type i = 0; i < nc; i++)
     for (octave_idx_type j = 0; j < nc; j++)
       retval(i,j) = tmpMat(invpvec(i),invpvec(j));
 
-  OCTAVE_QUIT;
-
-  for (octave_idx_type i = 0; i < nc; i++)
-    iperm(i) = i;  // identity permutation
-
-  // trailing permutations must be done in reverse order
-  for (octave_idx_type i = nc - 1; i >= ihi; i--)
-    {
-      octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
-      octave_idx_type tmp = iperm(i);
-      iperm(i) = iperm(swapidx);
-      iperm(swapidx) = tmp;
-    }
-
-  // construct inverse balancing permutation vector
-  for (octave_idx_type i = 0; i < nc; i++)
-    invpvec(iperm(i)) = i;     // Thanks to R. A. Lippert for this method
-
-  OCTAVE_QUIT;
- 
-  tmpMat = retval;
-  for (octave_idx_type i = 0; i < nc; i++)
-    for (octave_idx_type j = 0; j < nc; j++)
-      retval(i,j) = tmpMat(invpvec(i),invpvec(j));
-
-  // Reverse preconditioning step 1: fix trace normalization.
-  
+  // Reverse preconditioning step 1: fix trace normalization. 
   if (trshift > 0.0)
     retval = exp (trshift) * retval;