Mercurial > hg > octave-nkf

changeset 11842:0b9c56b6bf0e release-3-0-x

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
partially sync Matrix::expm and ComplexMatrix::expm with development repo
author Jaroslav Hajek <highegg@gmail.com>
date Fri, 19 Sep 2008 11:29:51 +0200
parents 43fccbab412a
children d1b8260dbc76
files liboctave/CMatrix.cc liboctave/dMatrix.cc
diffstat 2 files changed, 22 insertions(+), 10 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -2896,6 +2896,9 @@
 
   if (sqpow > 0)
     {
+      if (sqpow > 1023)
+	sqpow = 1023;
+
       double scale_factor = 1.0;
       for (octave_idx_type i = 0; i < sqpow; i++)
 	scale_factor *= 2.0;
@@ -2942,8 +2945,11 @@
 
   // Compute pade approximation = inverse (dpp) * npp.
 
-  retval = dpp.solve (npp);
-	
+  retval = dpp.solve (npp, info);
+
+  if (info < 0)
+    return retval;
+
   // Reverse preconditioning step 3: repeated squaring.
 
   while (sqpow)
@@ -2977,12 +2983,13 @@
     }
 
   // construct inverse balancing permutation vector
+  Array<octave_idx_type> 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));
@@ -3002,13 +3009,12 @@
     }
 
   // construct inverse balancing permutation vector
-  Array<octave_idx_type> 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;
+  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));
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -2523,10 +2523,13 @@
   
   if (sqpow > 0)
     {
+      if (sqpow > 1023)
+	sqpow = 1023;
+
       double scale_factor = 1.0;
       for (octave_idx_type i = 0; i < sqpow; i++)
 	scale_factor *= 2.0;
-  
+
       m = m / scale_factor;
     }
   
@@ -2570,7 +2573,10 @@
   // Compute pade approximation = inverse (dpp) * npp.
 
   retval = dpp.solve (npp, info);
-  
+
+  if (info < 0)
+    return retval;
+
   // Reverse preconditioning step 3: repeated squaring.
   
   while (sqpow)
@@ -2602,12 +2608,13 @@
     }
 
   // construct inverse balancing permutation vector
+  Array<octave_idx_type> 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;
+  Matrix 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));
@@ -2627,13 +2634,12 @@
     }
 
   // construct inverse balancing permutation vector
-  Array<octave_idx_type> 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;
  
-  Matrix tmpMat = retval;
+  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));