# HG changeset patch
# User Jaroslav Hajek
# Date 1221643678 -7200
# Node ID 87c7705d0be8c4268bddd6e51c24f84dc3721aef
# Parent fc45357bf50ccdb48fceb36ef92db2c5e19dc671
simplify expm inverse permutations
diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -3036,8 +3036,8 @@
for (octave_idx_type i = 0; i < nc; i++)
iperm(i) = i; // initialize to identity permutation
- // leading permutations in forward order
- for (octave_idx_type i = 0; i < (ilo-1); i++)
+ // trailing permutations must be done in reverse order
+ for (octave_idx_type i = nc - 1; i >= ihi; i--)
{
octave_idx_type swapidx = static_cast (dpermute(i)) - 1;
octave_idx_type tmp = iperm(i);
@@ -3045,25 +3045,8 @@
iperm(swapidx) = tmp;
}
- // construct inverse balancing permutation vector
- Array 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--)
+ // leading permutations in forward order
+ for (octave_idx_type i = 0; i < (ilo-1); i++)
{
octave_idx_type swapidx = static_cast (dpermute(i)) - 1;
octave_idx_type tmp = iperm(i);
@@ -3072,17 +3055,19 @@
}
// construct inverse balancing permutation vector
+ Array 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;
- tmpMat = retval;
+ 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));
// 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
@@ -0,0 +1,5 @@
+2008-09-12 Marco Caliari
+
+ * dMatrix.cc (Matrix::expm): Simplify reverse permutation.
+ * CMatrix.cc (ComplexMatrix::expm): Likewise.
+
diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -2668,12 +2668,21 @@
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 (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 (dpermute(i)) - 1;
octave_idx_type tmp = iperm(i);
- iperm(i) = iperm (swapidx);
+ iperm(i) = iperm(swapidx);
iperm(swapidx) = tmp;
}
@@ -2689,32 +2698,8 @@
for (octave_idx_type j = 0; j < nc; j++)
retval(i,j) = tmpMat(invpvec(i),invpvec(j));
- OCTAVE_QUIT;
+ // Reverse preconditioning step 1: fix trace normalization.
- 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 (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.
if (trshift > 0.0)
retval = exp (trshift) * retval;