# HG changeset patch
# User jwe
# Date 948492793 0
# Node ID a2dc6de198f9840183d2657e0e9be45bfc200502
# Parent  63378e4a34f23fa1fd2fa47983afe2914af36fc7
[project @ 2000-01-21 22:13:13 by jwe]

diff -r 63378e4a34f2 -r a2dc6de198f9 liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc	Fri Jan 21 19:31:51 2000 +0000
+++ b/liboctave/CMatrix.cc	Fri Jan 21 22:13:13 2000 +0000
@@ -1606,11 +1606,13 @@
   Complex *mp = m.fortran_vec ();
 
   int info, ilo, ihi,ilos,ihis;
-  Array<double> dpermute(nc);
-  Array<double> dscale(nc);
-
-  // FIXME: should pass job as a parameter in expm
-  char job = 'P';       // Permute first
+  Array<double> dpermute (nc);
+  Array<double> dscale (nc);
+
+  // XXX FIXME XXX -- should pass job as a parameter in expm
+
+  // Permute first
+  char job = 'P';
   F77_XFCN (zgebal, ZGEBAL, (&job, nc, mp, nc, ilo, ihi,
             dpermute.fortran_vec(), info, 1L, 1L));
 
@@ -1620,7 +1622,8 @@
       return retval;
     }
 
-  job = 'S';            // then scale
+  // then scale
+  job = 'S';
   F77_XFCN (zgebal, ZGEBAL, (&job, nc, mp, nc, ilos, ihis,
             dscale.fortran_vec(), info, 1L, 1L));
 
@@ -1701,43 +1704,42 @@
   // Done in two steps: inverse scaling, then inverse permutation
 
   // inverse scaling (diagonal transformation)
-  int ii, jj;
-  for(ii = 0; ii < nc ; ii++)
-    for(jj=0; jj < nc ; jj++)
-       retval(ii,jj) *= dscale(ii)/dscale(jj);
+  for (int i = 0; i < nc; i++)
+    for (int j = 0; j < nc; j++)
+       retval(i,j) *= dscale(i) / dscale(j);
 
   // construct balancing permutation vector
-  Array<int> ipermute(nc);
-  for(ii=0 ; ii < nc ; ii++)
-    ipermute(ii) = ii;  // initialize to identity permutation
+  Array<int> ipermute (nc);
+  for (int i = 0; i < nc; i++)
+    ipermute(i) = i;  // initialize to identity permutation
 
   // leading permutations in forward order
-  for( ii = 0 ; ii < (ilo-1) ; ii++)
-  {
-    int swapidx = ( (int) dpermute(ii) ) -1;
-    int tmp = ipermute(ii);
-    ipermute(ii) = ipermute( swapidx );
-    ipermute(swapidx) = tmp;
-  }
+  for (int i = 0; i < (ilo-1); i++)
+    {
+      int swapidx = static_cast<int> (dpermute(i)) - 1;
+      int tmp = ipermute(i);
+      ipermute(i) = ipermute(swapidx);
+      ipermute(swapidx) = tmp;
+    }
 
   // trailing permutations must be done in reverse order
-  for( ii = nc-1 ; ii >= ihi ; ii--)
-  {
-    int swapidx = ( (int) dpermute(ii) ) -1;
-    int tmp = ipermute(ii);
-    ipermute(ii) = ipermute( swapidx );
-    ipermute(swapidx) = tmp;
-  }
+  for (int i = nc - 1; i >= ihi; i--)
+    {
+      int swapidx = static_cast<int> (dpermute(i)) - 1;
+      int tmp = ipermute(i);
+      ipermute(i) = ipermute(swapidx);
+      ipermute(swapidx) = tmp;
+    }
 
   // construct inverse balancing permutation vector
-  Array<int> invpvec(nc);
-  for( ii = 0 ; ii < nc ; ii++)
-    invpvec(ipermute(ii)) = ii ;     // Thanks to R. A. Lippert for this method
+  Array<int> invpvec (nc);
+  for (int i = 0; i < nc; i++)
+    invpvec(ipermute(i)) = i;     // Thanks to R. A. Lippert for this method
 
   ComplexMatrix tmpMat = retval;
-  for( ii = 0 ; ii < nc ; ii ++)
-    for( jj= 0 ; jj < nc ; jj++ )
-      retval(ii,jj) = tmpMat(invpvec(ii),invpvec(jj));
+  for (int i = 0; i < nc; i++)
+    for (int j = 0; j < nc; j++)
+      retval(i,j) = tmpMat(invpvec(i),invpvec(j));
 
   // Reverse preconditioning step 1: fix trace normalization.
 
diff -r 63378e4a34f2 -r a2dc6de198f9 liboctave/ChangeLog
--- a/liboctave/ChangeLog	Fri Jan 21 19:31:51 2000 +0000
+++ b/liboctave/ChangeLog	Fri Jan 21 22:13:13 2000 +0000
@@ -1,3 +1,10 @@
+2000-01-21  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* CMatrix.cc (ComplexMatrix::expm): Apply permutation and scaling
+	operations directly in step 2 and reverse step 2.
+	* dMatrix.cc (Matrix::expm): Apply permutation and scaling
+	operations directly in step 2 and reverse step 2.
+
 2000-01-20  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* oct-time.h, oct-time.cc (octave_strptime): New class.
diff -r 63378e4a34f2 -r a2dc6de198f9 liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc	Fri Jan 21 19:31:51 2000 +0000
+++ b/liboctave/dMatrix.cc	Fri Jan 21 22:13:13 2000 +0000
@@ -1356,16 +1356,18 @@
 
   double *p_m = m.fortran_vec ();
 
-  int info, ilo, ihi,ilos,ihis, ii, jj;
-  Array<double> dpermute(nc);
-  Array<double> dscale(nc);
+  int info, ilo, ihi, ilos, ihis;
+  Array<double> dpermute (nc);
+  Array<double> dscale (nc);
   double *dp;
 
-  char job = 'P';       // permutation first
+  // permutation first
+  char job = 'P';
   dp = dpermute.fortran_vec();
   F77_XFCN (dgebal, DGEBAL, (&job, nc, p_m, nc, ilo, ihi, dp, info, 1L, 1L));
 
-  job = 'S';       // then scaling
+  // then scaling
+  job = 'S';
   dp = dscale.fortran_vec();
   F77_XFCN (dgebal, DGEBAL, (&job, nc, p_m, nc, ilos, ihis, dp, info, 1L, 1L));
 
@@ -1445,41 +1447,42 @@
   
   // Reverse preconditioning step 2: inverse balancing.
   // apply inverse scaling to computed exponential
-  for(ii = 0; ii < nc ; ii++)
-    for(jj=0; jj < nc ; jj++)
-       retval(ii,jj) *= dscale(ii)/dscale(jj);
+  for (int i = 0; i < nc; i++)
+    for (int j = 0; j < nc; j++)
+       retval(i,j) *= dscale(i) / dscale(j);
 
   // construct balancing permutation vector
-  Array<int> ipermute(nc);
-  for(ii=0 ; ii < nc ; ii++) ipermute(ii) = ii;  // identity permutation
+  Array<int> ipermute (nc);
+  for (int i = 0; i < nc; i++)
+    ipermute(i) = i;  // identity permutation
 
   // leading permutations in forward order
-  for( ii = 0 ; ii < (ilo-1) ; ii++)
-  {
-    int swapidx = ( (int) dpermute(ii) ) -1;
-    int tmp = ipermute(ii);
-    ipermute(ii) = ipermute( swapidx );
-    ipermute(swapidx) = tmp;
-  }
+  for (int i = 0; i < (ilo-1); i++)
+    {
+      int swapidx = static_cast<int> (dpermute(i)) - 1;
+      int tmp = ipermute(i);
+      ipermute(i) = ipermute (swapidx);
+      ipermute(swapidx) = tmp;
+    }
 
   // trailing permutations must be done in reverse order
-  for( ii = nc-1 ; ii >= ihi ; ii--)
-  {
-    int swapidx = ( (int) dpermute(ii) ) -1;
-    int tmp = ipermute(ii);
-    ipermute(ii) = ipermute( swapidx );
-    ipermute(swapidx) = tmp;
-  }
+  for (int i = nc - 1; i >= ihi; i--)
+    {
+      int swapidx = static_cast<int> (dpermute(i)) - 1;
+      int tmp = ipermute(i);
+      ipermute(i) = ipermute(swapidx);
+      ipermute(swapidx) = tmp;
+    }
 
   // construct inverse balancing permutation vector
-  Array<int> invpvec(nc);
-  for( ii = 0 ; ii < nc ; ii++)
-    invpvec(ipermute(ii)) = ii ;     // Thanks to R. A. Lippert for this method
+  Array<int> invpvec (nc);
+  for (int i = 0; i < nc; i++)
+    invpvec(ipermute(i)) = i;     // Thanks to R. A. Lippert for this method
  
   Matrix tmpMat = retval;
-  for( ii = 0 ; ii < nc ; ii ++)
-    for( jj= 0 ; jj < nc ; jj++ )
-      retval(ii,jj) = tmpMat(invpvec(ii),invpvec(jj));
+  for (int i = 0; i < nc; i++)
+    for (int j = 0; j < nc; j++)
+      retval(i,j) = tmpMat(invpvec(i),invpvec(j));
 
   // Reverse preconditioning step 1: fix trace normalization.