# HG changeset patch
# User jwe
# Date 1201163496 0
# Node ID a7a987b229b7b0e16b2ae8748e147c2caffa815d
# Parent  2da391ef7f55f26abbc157977a96ccf1040f026f
[project @ 2008-01-24 08:31:36 by jwe]

diff -r 2da391ef7f55 -r a7a987b229b7 liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc	Thu Jan 24 07:40:58 2008 +0000
+++ b/liboctave/CMatrix.cc	Thu Jan 24 08:31:36 2008 +0000
@@ -2938,6 +2938,23 @@
       iperm(swapidx) = tmp;
     }
 
+  // construct inverse balancing permutation vector
+  Array<octave_idx_type> invpvec (nc);
+  for (octave_idx_type i = 0; i < nc; i++)
+    invpvec(iperm(i)) = i;     // Thanks to R. A. Lippert for this method
+
+  OCTAVE_QUIT;
+
+  ComplexMatrix 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;  // initialize to identity permutation
+
   // trailing permutations must be done in reverse order
   for (octave_idx_type i = nc - 1; i >= ihi; i--)
     {
@@ -2948,13 +2965,12 @@
     }
 
   // construct inverse balancing permutation vector
-  Array<octave_idx_type> invpvec (nc);
   for (octave_idx_type i = 0; i < nc; i++)
     invpvec(iperm(i)) = i;     // Thanks to R. A. Lippert for this method
 
   OCTAVE_QUIT;
 
-  ComplexMatrix tmpMat = retval;
+  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));
diff -r 2da391ef7f55 -r a7a987b229b7 liboctave/ChangeLog
--- a/liboctave/ChangeLog	Thu Jan 24 07:40:58 2008 +0000
+++ b/liboctave/ChangeLog	Thu Jan 24 08:31:36 2008 +0000
@@ -1,3 +1,8 @@
+2008-01-18  Marco Caliari  <marco.caliari@univr.it>
+
+	* dMatrix.cc (Matrix::expm): Correctly perform reverse permutation.
+	* CMatrix.cc (ComplexMatrix::expm): Likewise.
+
 2008-01-22  Michael Goffioul  <michael.goffioul@gmail.com>
 
 	* oct-time.cc (octave_base_tim::init): Validate pointer argument;
diff -r 2da391ef7f55 -r a7a987b229b7 liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc	Thu Jan 24 07:40:58 2008 +0000
+++ b/liboctave/dMatrix.cc	Thu Jan 24 08:31:36 2008 +0000
@@ -2556,6 +2556,23 @@
       iperm(swapidx) = tmp;
     }
 
+  // construct inverse balancing permutation vector
+  Array<octave_idx_type> invpvec (nc);
+  for (octave_idx_type i = 0; i < nc; i++)
+    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--)
     {
@@ -2566,13 +2583,12 @@
     }
 
   // construct inverse balancing permutation vector
-  Array<octave_idx_type> invpvec (nc);
   for (octave_idx_type i = 0; i < nc; i++)
     invpvec(iperm(i)) = i;     // Thanks to R. A. Lippert for this method
 
   OCTAVE_QUIT;
  
-  Matrix tmpMat = retval;
+  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));