Expm small patches

[Marco Caliari]

Hi.

I have two small patches for expm in liboctave/dMatrix.cc and 
liboctave/CMatrix.cc. They essentially substitute elem with fortran_vec 
everywhere.

Best regards,

Marco
-------------- next part --------------
--- liboctave/dMatrix.cc.orig	2007-12-14 10:06:10.000000000 +0100
+++ liboctave/dMatrix.cc	2007-12-14 10:07:33.000000000 +0100
@@ -2384,9 +2384,7 @@
   Matrix retval;
 
   Matrix m = *this;
-
-  if (numel () == 1)
-    return Matrix (1, 1, exp (m(0)));
+  double *mp = m.fortran_vec ();
 
   octave_idx_type nc = columns ();
 
@@ -2397,21 +2395,25 @@
   volatile double trshift = 0.0;
 
   for (octave_idx_type i = 0; i < nc; i++)
-    trshift += m.elem (i, i);
+    {
+      octave_idx_type k = i * nc + i;
+      trshift += mp [k];
+    }
 
   trshift /= nc;
 
   if (trshift > 0.0)
     {
       for (octave_idx_type i = 0; i < nc; i++)
-	m.elem (i, i) -= trshift;
+	{
+	  octave_idx_type k = i * nc + i;
+	  mp [k] -= trshift;
+	}
     }
 
   // Preconditioning step 2: balancing; code follows development
   // in AEPBAL
 
-  double *p_m = m.fortran_vec ();
-
   octave_idx_type info, ilo, ihi, ilos, ihis;
   Array<double> dpermute (nc);
   Array<double> dscale (nc);
@@ -2419,14 +2421,14 @@
   // permutation first
   char job = 'P';
   F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
-			     nc, p_m, nc, ilo, ihi,
+			     nc, mp, nc, ilo, ihi,
 			     dpermute.fortran_vec (), info
 			     F77_CHAR_ARG_LEN (1)));
 
   // then scaling
   job = 'S';
   F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
-			     nc, p_m, nc, ilos, ihis,
+			     nc, mp, nc, ilos, ihis,
 			     dscale.fortran_vec (), info
 			     F77_CHAR_ARG_LEN (1)));
 
@@ -2479,32 +2481,31 @@
   
   // Now powers a^8 ... a^1.
   
-  octave_idx_type minus_one_j = -1;
+  octave_idx_type minus_one_j = 1;
   for (octave_idx_type j = 7; j >= 0; j--)
     {
       for (octave_idx_type i = 0; i < nc; i++)
 	{
 	  octave_idx_type k = i * nc + i;
-	  pnpp[k] += padec[j];
-	  pdpp[k] += minus_one_j * padec[j];
+	  pnpp [k] += padec [j];
+	  pdpp [k] += minus_one_j * padec [j];
 	}      
-
       npp = m * npp;
       pnpp = npp.fortran_vec ();
 
       dpp = m * dpp;
       pdpp = dpp.fortran_vec ();
 
-      minus_one_j *= -1;
+      minus_one_j = -minus_one_j;
     }
   
   // Zero power.
   
-  dpp = -dpp;
   for (octave_idx_type j = 0; j < nc; j++)
     {
-      npp.elem (j, j) += 1.0;
-      dpp.elem (j, j) += 1.0;
+      octave_idx_type k = j * nc + j;
+      pnpp [k] += 1.0;
+      pdpp [k] += 1.0;
     }
   
   // Compute pade approximation = inverse (dpp) * npp.
-------------- next part --------------
--- liboctave/CMatrix.cc.orig	2007-12-14 10:06:18.000000000 +0100
+++ liboctave/CMatrix.cc	2007-12-14 10:08:06.000000000 +0100
@@ -2756,6 +2756,7 @@
   ComplexMatrix retval;
 
   ComplexMatrix m = *this;
+  Complex *mp = m.fortran_vec ();
 
   octave_idx_type nc = columns ();
 
@@ -2766,8 +2767,10 @@
   Complex trshift = 0.0;
 
   for (octave_idx_type i = 0; i < nc; i++)
-    trshift += m.elem (i, i);
-
+    {
+      octave_idx_type k = i * nc + i;
+      trshift += mp [k];
+    }
   trshift /= nc;
 
   if (trshift.real () < 0.0)
@@ -2778,13 +2781,14 @@
     }
 
   for (octave_idx_type i = 0; i < nc; i++)
-    m.elem (i, i) -= trshift;
+    {
+      octave_idx_type k = i * nc + i;
+      mp [k] -= trshift;
+    }
 
   // Preconditioning step 2: eigenvalue balancing.
   // code follows development in AEPBAL
 
-  Complex *mp = m.fortran_vec ();
-
   octave_idx_type info, ilo, ihi,ilos,ihis;
   Array<double> dpermute (nc);
   Array<double> dscale (nc);
@@ -2833,8 +2837,9 @@
       return retval;
     }
 
-  int sqpow = (inf_norm > 0.0
-	       ? static_cast<int> (1.0 + log (inf_norm) / log (2.0)) : 0);
+  octave_idx_type sqpow = static_cast<octave_idx_type> (inf_norm > 0.0
+		     ? (1.0 + log (inf_norm) / log (2.0))
+		     : 0.0);
 
   // Check whether we need to square at all.
 
@@ -2859,32 +2864,31 @@
 
   // Now powers a^8 ... a^1.
 
-  int minus_one_j = -1;
+  octave_idx_type minus_one_j = 1;
   for (octave_idx_type j = 7; j >= 0; j--)
     {
       for (octave_idx_type i = 0; i < nc; i++)
 	{
 	  octave_idx_type k = i * nc + i;
-	  pnpp[k] += padec[j];
-	  pdpp[k] += minus_one_j * padec[j];
+	  pnpp [k] += padec [j];
+	  pdpp [k] += minus_one_j * padec [j];
 	}      
-
       npp = m * npp;
       pnpp = npp.fortran_vec ();
 
       dpp = m * dpp;
       pdpp = dpp.fortran_vec ();
 
-      minus_one_j *= -1;
+      minus_one_j = -minus_one_j;
     }
 
   // Zero power.
 
-  dpp = -dpp;
   for (octave_idx_type j = 0; j < nc; j++)
     {
-      npp.elem (j, j) += 1.0;
-      dpp.elem (j, j) += 1.0;
+      octave_idx_type k = j * nc + j;
+      pnpp [k] += 1.0;
+      pdpp [k] += 1.0;
     }
 
   // Compute pade approximation = inverse (dpp) * npp.