Mercurial > hg > octave-nkf

changeset 11841:43fccbab412a release-3-0-x

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2008-01-24 08:31:36 by jwe]
author jwe
date Wed, 17 Sep 2008 11:06:18 +0200
parents b160651f8a21
children 0b9c56b6bf0e
files liboctave/CMatrix.cc liboctave/ChangeLog liboctave/dMatrix.cc
diffstat 3 files changed, 37 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -2976,6 +2976,22 @@
       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));
+
+  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--)
     {
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -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-09-12  Jaroslav Hajek  <highegg@gmail.com>
 
 	* oct-inttypes.h (pow (const octave_int<T>&, const octave_int<T>&)): 
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -2601,6 +2601,22 @@
       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));
+
+  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--)
     {