# HG changeset patch
# User jwe
# Date 1200644724 0
# Node ID f9df7f7520e7c981c8b4712ca16958ca7dd0f80b
# Parent  c1a3d6c7d2fbd8f152fa8b6e7814f4835d82cd35
[project @ 2008-01-18 08:25:24 by jwe]

diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -2750,6 +2750,14 @@
   1.9270852604185938e-9,
 };
 
+static void
+solve_singularity_warning (double rcond)
+{
+  (*current_liboctave_warning_handler) 
+    ("singular matrix encountered in expm calculation, rcond = %g",
+     rcond);
+}
+
 ComplexMatrix
 ComplexMatrix::expm (void) const
 {
@@ -2843,6 +2851,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;
@@ -2889,8 +2900,12 @@
 
   // Compute pade approximation = inverse (dpp) * npp.
 
-  retval = dpp.solve (npp);
-	
+  double rcond;
+  retval = dpp.solve (npp, info, rcond, solve_singularity_warning);
+
+  if (info < 0)
+    return retval;
+
   // Reverse preconditioning step 3: repeated squaring.
 
   while (sqpow)
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,15 @@
+2008-01-18  John W. Eaton  <jwe@octave.org>
+
+	* dMatrix.cc (solve_singularity_warning): New function.
+	(Matrix::expm): Pass pointer to solve_singularity_warning to
+	Matrix::solve method.  Exit early if Matrix::solve fails.
+	Limit sqpow value to avoid overflowing scale factor.
+	* CMatrix.cc (solve_singularity_warning): New function.
+	(ComplexMatrix::expm): Pass pointer to solve_singularity_warning to
+	ComplexMatrix::solve method.  Exit early if ComplexMatrix::solve fails.
+	Limit sqpow value to avoid overflowing scale factor.
+	From Marco Caliari <mcaliari@math.unipd.it>.
+
 2008-01-10  Kim Hansen  <kimhanse@gmail.com>
 
 	* Sparse.cc: New tests for slicing of sparse matrices.
diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -2378,6 +2378,14 @@
   1.9270852604185938e-9,
 };
 
+static void
+solve_singularity_warning (double rcond)
+{
+  (*current_liboctave_warning_handler) 
+    ("singular matrix encountered in expm calculation, rcond = %g",
+     rcond);
+}
+
 Matrix
 Matrix::expm (void) const
 {
@@ -2463,10 +2471,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;
     }
   
@@ -2509,8 +2520,12 @@
   
   // Compute pade approximation = inverse (dpp) * npp.
 
-  retval = dpp.solve (npp, info);
-  
+  double rcond;
+  retval = dpp.solve (npp, info, rcond, solve_singularity_warning);
+
+  if (info < 0)
+    return retval;
+
   // Reverse preconditioning step 3: repeated squaring.
   
   while (sqpow)