--- 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.