Add exception handler to memory allocation

11 files changed, 99 insertions, 94 deletions

diff --git a/Eigen/src/OrderingMethods/Eigen_Colamd.h b/Eigen/src/OrderingMethods/Eigen_Colamd.h
index 39701d0af..0af137d54 100644
--- a/Eigen/src/OrderingMethods/Eigen_Colamd.h
+++ b/Eigen/src/OrderingMethods/Eigen_Colamd.h
@@ -155,7 +155,6 @@
 #endif /* MATLAB_MEX_FILE */
 
     // == Row and Column structures ==
-
 typedef struct EIGEN_Colamd_Col_struct
 {
     int start ;   /* index for A of first row in this column, or EIGEN_DEAD */
@@ -248,11 +247,9 @@ void eigen_colamd_report (int stats [EIGEN_COLAMD_STATS]);
 
 int eigen_init_rows_cols (int n_row, int n_col, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col col [], int A [], int p [], int stats[EIGEN_COLAMD_STATS] ); 
 
-void eigen_init_scoring (int n_row, int n_col, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col Col [], int A [], int head [], 
-                    double knobs[EIGEN_COLAMD_KNOBS], int *p_n_row2, int *p_n_col2, int *p_max_deg);
+void eigen_init_scoring (int n_row, int n_col, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col Col [], int A [], int head [], double knobs[EIGEN_COLAMD_KNOBS], int *p_n_row2, int *p_n_col2, int *p_max_deg);
 
-int eigen_find_ordering (int n_row, int n_col, int Alen, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col Col [], int A [], int head [], 
-                    int n_col2, int max_deg, int pfree);
+int eigen_find_ordering (int n_row, int n_col, int Alen, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col Col [], int A [], int head [], int n_col2, int max_deg, int pfree);
 
 void eigen_order_children (int n_col, EIGEN_Colamd_Col Col [], int p []);
 
@@ -2514,5 +2511,4 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E
 }
 
 #endif /* NDEBUG */
-
 #endif
diff --git a/Eigen/src/OrderingMethods/Ordering.h b/Eigen/src/OrderingMethods/Ordering.h
index cbd2e5d34..47cd6f169 100644
--- a/Eigen/src/OrderingMethods/Ordering.h
+++ b/Eigen/src/OrderingMethods/Ordering.h
@@ -124,9 +124,6 @@ class COLAMDOrdering
     typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; 
     typedef Matrix<Index, Dynamic, 1> IndexVector; 
     /** Compute the permutation vector form a sparse matrix */
-    
-
-
     template <typename MatrixType>
     void operator() (const MatrixType& mat, PermutationType& perm)
     {
@@ -152,9 +149,6 @@ class COLAMDOrdering
         for (int i = 0; i < n; i++) perm.indices()(p(i)) = i;
         
     }
-    
-  private:
-    
  
 };
 
diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
index bb1decc4c..25fad0f29 100644
--- a/Eigen/src/SparseLU/SparseLU.h
+++ b/Eigen/src/SparseLU/SparseLU.h
@@ -339,9 +339,6 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
   
   //TODO  It is possible as in SuperLU to compute row and columns scaling vectors to equilibrate the matrix mat.
   
-  // Compute the fill-reducing ordering  
-  // TODO Currently, the only  available ordering method is AMD. 
-  
   OrderingType ord; 
   ord(mat,m_perm_c);
   //FIXME Check the right semantic behind m_perm_c
diff --git a/Eigen/src/SparseLU/SparseLU_Coletree.h b/Eigen/src/SparseLU/SparseLU_Coletree.h
index 1329d383f..142f4995e 100644
--- a/Eigen/src/SparseLU/SparseLU_Coletree.h
+++ b/Eigen/src/SparseLU/SparseLU_Coletree.h
@@ -94,7 +94,6 @@ int LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
   int rset, cset, rroot; 
   for (col = 0; col < nc; col++) 
   {
-//     cset = pp(col) = col; // Initially, each element is in its own set //FIXME
     pp(col) = col; 
     cset = col; 
     root(cset) = col; 
@@ -108,7 +107,6 @@ int LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
       if (rroot != col) 
       {
         parent(rroot) = col; 
-//         cset = pp(cset) = rset; // Get the union of cset and rset  //FIXME
         pp(cset) = rset; 
         cset = rset; 
         root(cset) = col; 
diff --git a/Eigen/src/SparseLU/SparseLU_Matrix.h b/Eigen/src/SparseLU/SparseLU_Matrix.h
index 90a0f2740..9f2dcaa56 100644
--- a/Eigen/src/SparseLU/SparseLU_Matrix.h
+++ b/Eigen/src/SparseLU/SparseLU_Matrix.h
@@ -192,7 +192,6 @@ class SuperNodalMatrix
   protected:
     Index m_row; // Number of rows
     Index m_col; // Number of columns 
-//     Index m_nnz; // Number of nonzero values 
     Index m_nsuper; // Number of supernodes 
     Scalar* m_nzval; //array of nonzero values packed by column
     Index* m_nzval_colptr; //nzval_colptr[j] Stores the location in nzval[] which starts column j 
diff --git a/Eigen/src/SparseLU/SparseLU_Memory.h b/Eigen/src/SparseLU/SparseLU_Memory.h
index a17079199..7a2ab93df 100644
--- a/Eigen/src/SparseLU/SparseLU_Memory.h
+++ b/Eigen/src/SparseLU/SparseLU_Memory.h
@@ -78,41 +78,82 @@ int  expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_ex
   
   VectorType old_vec; // Temporary vector to hold the previous values   
   if (nbElts > 0 )
-    old_vec = vec.segment(0,nbElts); // old_vec should be of size nbElts... to be checked
+    old_vec = vec.segment(0,nbElts); 
   
-  //expand the current vector //FIXME Should be in a try ... catch region
-  vec.resize(new_len); 
-  /* 
-   * Test if the memory has been well allocated  
-   * otherwise reduce the size and try to reallocate
-   * copy data from previous vector (if exists) to the newly allocated vector 
-   */
-  if ( num_expansions != 0 ) // The memory has been expanded before
+  //Allocate or expand the current vector
+  try 
+  {
+    vec.resize(new_len); 
+  }
+  catch(std::bad_alloc& )
   {
-    int tries = 0; 
-    if (keep_prev) 
+    if ( !num_expansions )
     {
-      if (!vec.size()) return new_len ;
+      // First time to allocate from LUMemInit()
+      throw; // Pass the exception to LUMemInit() which has a try... catch block
+    }
+    if (keep_prev)
+    {
+      // In this case, the memory length should not not be reduced
+      return new_len;
     }
     else 
     {
-      while (!vec.size())
+      // Reduce the size and increase again 
+      int tries = 0; // Number of attempts
+      do 
       {
-       // Reduce the size and allocate again
-        if ( ++tries > 10) return new_len;  
-        alpha = LU_Reduce(alpha); 
+        alpha = LU_Reduce(alpha);
         new_len = alpha * length ; 
-        vec.resize(new_len); //FIXME Should be in a try catch section
-      }
-    } // end allocation 
-    
-    //Copy the previous values to the newly allocated space 
-    if (nbElts > 0)
-      vec.segment(0, nbElts) = old_vec;   
-  } // end expansion 
+        try
+        {
+          vec.resize(new_len); 
+        }
+        catch(std::bad_alloc& )
+        {
+          tries += 1; 
+          if ( tries > 10) return new_len; 
+        }
+      } while (!vec.size());
+    }
+  }
+  //Copy the previous values to the newly allocated space 
+  if (nbElts > 0)
+    vec.segment(0, nbElts) = old_vec;   
+   
+  
   length  = new_len;
   if(num_expansions) ++num_expansions;
   return 0; 
+  
+  /* 
+   * Test if the memory has been well allocated  
+   * otherwise reduce the size and try to reallocate
+   * copy data from previous vector (if exists) to the newly allocated vector 
+   */
+//   if ( num_expansions != 0 ) // The memory has been expanded before
+//   {
+//     int tries = 0; 
+//     if (keep_prev) 
+//     {
+//       if (!vec.size()) return new_len ;
+//     }
+//     else 
+//     {
+//       while (!vec.size())
+//       {
+//        // Reduce the size and allocate again
+//         if ( ++tries > 10) return new_len;  
+//         alpha = LU_Reduce(alpha); 
+//         new_len = alpha * length ; 
+//         vec.resize(new_len); //FIXME Should be in a try catch section
+//       }
+//     } // end allocation 
+//     
+//     //Copy the previous values to the newly allocated space 
+//     if (nbElts > 0)
+//       vec.segment(0, nbElts) = old_vec;   
+//   } // end expansion
 }
 
 /**
@@ -122,8 +163,8 @@ int  expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_ex
  * \param annz number of initial nonzeros in the matrix 
  * \param lwork  if lwork=-1, this routine returns an estimated size of the required memory
  * \param glu persistent data to facilitate multiple factors : will be deleted later ??
- * \return an estimated size of the required memory if lwork = -1; otherwise, return the size of actually allocated when memory allocation failed 
- * NOTE Unlike SuperLU, this routine does not support successive factorization with the same pattern and the row permutation
+ * \return an estimated size of the required memory if lwork = -1; otherwise, return the size of actually allocated memory when allocation failed, and 0 on success
+ * NOTE Unlike SuperLU, this routine does not support successive factorization with the same pattern and the same row permutation
  */
 template <typename IndexVector,typename ScalarVector>
 int LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size,  LU_GlobalLU_t<IndexVector,ScalarVector>& glu)
@@ -159,27 +200,26 @@ int LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size,
   glu.xusub.resize(n+1);
 
   // Reserve memory for L/U factors
-  expand<ScalarVector>(glu.lusup, nzlumax, 0, 0, num_expansions); 
-  expand<ScalarVector>(glu.ucol,nzumax, 0, 0, num_expansions); 
-  expand<IndexVector>(glu.lsub,nzlmax, 0, 0, num_expansions); 
-  expand<IndexVector>(glu.usub,nzumax, 0, 1, num_expansions);
-  
-  // Check if the memory is correctly allocated, 
-  // FIXME Should be a try... catch section here 
-  while ( !glu.lusup.size() || !glu.ucol.size() || !glu.lsub.size() || !glu.usub.size())
+  do 
   {
-    //Reduce the estimated size and retry
-    nzlumax /= 2;
-    nzumax /= 2;
-    nzlmax /= 2;
-    
-    if (nzlumax < annz ) return nzlumax; 
+    try
+    {
+      expand<ScalarVector>(glu.lusup, nzlumax, 0, 0, num_expansions); 
+      expand<ScalarVector>(glu.ucol,nzumax, 0, 0, num_expansions); 
+      expand<IndexVector>(glu.lsub,nzlmax, 0, 0, num_expansions); 
+      expand<IndexVector>(glu.usub,nzumax, 0, 1, num_expansions); 
+    }
+    catch(std::bad_alloc& )
+    {
+      //Reduce the estimated size and retry
+      nzlumax /= 2;
+      nzumax /= 2;
+      nzlmax /= 2;
+      if (nzlumax < annz ) return nzlumax; 
+    }
     
-    expand<ScalarVector>(glu.lusup, nzlumax, 0, 0, num_expansions); 
-    expand<ScalarVector>(glu.ucol, nzumax, 0, 0, num_expansions); 
-    expand<IndexVector>(glu.lsub, nzlmax, 0, 0, num_expansions); 
-    expand<IndexVector>(glu.usub, nzumax, 0, 1, num_expansions); 
-  }
+  } while (!glu.lusup.size() || !glu.ucol.size() || !glu.lsub.size() || !glu.usub.size());
+
   
   
   ++num_expansions;
@@ -207,23 +247,6 @@ int LUMemXpand(VectorType& vec, int& maxlen, int nbElts, LU_MemType memtype, int
   if (failed_size)
     return failed_size; 
   
-  // The following code  is not really needed since maxlen is passed by reference 
-  // and correspond to the appropriate field in glu
-//   switch ( mem_type ) {
-//     case LUSUP:
-//       glu.nzlumax = maxlen;
-//       break; 
-//     case UCOL:
-//       glu.nzumax = maxlen; 
-//       break; 
-//     case LSUB:
-//       glu.nzlmax = maxlen; 
-//       break;
-//     case USUB:
-//       glu.nzumax = maxlen; 
-//       break; 
-//   }
-  
   return 0 ; 
   
 }
diff --git a/Eigen/src/SparseLU/SparseLU_column_bmod.h b/Eigen/src/SparseLU/SparseLU_column_bmod.h
index 3042eb5f8..00787721b 100644
--- a/Eigen/src/SparseLU/SparseLU_column_bmod.h
+++ b/Eigen/src/SparseLU/SparseLU_column_bmod.h
@@ -133,7 +133,6 @@ int LU_column_bmod(const int jcol, const int nseg, BlockScalarVector& dense, Sca
       // Dense triangular solve -- start effective triangle
       luptr += nsupr * no_zeros + no_zeros; 
       // Form Eigen matrix and vector 
-//       std::cout<< "jcol " << jcol << " rows " << segsize << std::endl;
       Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) );
       VectorBlock<ScalarVector> u(tempv, 0, segsize);
       
diff --git a/Eigen/src/SparseLU/SparseLU_panel_dfs.h b/Eigen/src/SparseLU/SparseLU_panel_dfs.h
index 908ee67ac..70ea0f51f 100644
--- a/Eigen/src/SparseLU/SparseLU_panel_dfs.h
+++ b/Eigen/src/SparseLU/SparseLU_panel_dfs.h
@@ -92,7 +92,6 @@ void LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, Index
   int xdfs, maxdfs, kpar;
   
   // Initialize pointers 
-//   IndexVector& marker1 = marker.block(m, m); 
   VectorBlock<IndexVector> marker1(marker, m, m); 
   nseg = 0; 
   IndexVector& xsup = glu.xsup; 
diff --git a/Eigen/src/SparseLU/SparseLU_snode_bmod.h b/Eigen/src/SparseLU/SparseLU_snode_bmod.h
index 44438d037..fc8042f52 100644
--- a/Eigen/src/SparseLU/SparseLU_snode_bmod.h
+++ b/Eigen/src/SparseLU/SparseLU_snode_bmod.h
@@ -80,7 +80,6 @@ int LU_snode_bmod (const int jcol, const int fsupc, ScalarVector& dense, LU_Glob
     
     // Update the trailing part of the column jcol U(jcol:jcol+nrow, jcol) using L(jcol:jcol+nrow, fsupc:jcol) and U(fsupc:jcol)
     new (&A) Map<Matrix<Scalar,Dynamic,Dynamic>,0,OuterStride<> > ( &(lusup.data()[luptr+nsupc]), nrow, nsupc, OuterStride<>(nsupr) ); 
-//     Map<Matrix<Scalar,Dynamic,1> > l(&(lusup.data()[ufirst+nsupc], nrow); 
     VectorBlock<ScalarVector> l(lusup, ufirst+nsupc, nrow); 
     l = l - A * u; 
   }
diff --git a/bench/spbench/CMakeLists.txt b/bench/spbench/CMakeLists.txt
index 4b3c6f8e3..a093cc5d9 100644
--- a/bench/spbench/CMakeLists.txt
+++ b/bench/spbench/CMakeLists.txt
@@ -67,4 +67,4 @@ add_executable(spsolver sp_solver.cpp)
 target_link_libraries (spsolver ${SPARSE_LIBS})
 
 add_executable(test_sparseLU test_sparseLU.cpp)
-target_link_libraries (test_sparseLU ${SPARSE_LIBS})
\ No newline at end of file
+target_link_libraries (test_sparseLU ${SPARSE_LIBS})
diff --git a/bench/spbench/test_sparseLU.cpp b/bench/spbench/test_sparseLU.cpp
index 6fbf03454..31273add5 100644
--- a/bench/spbench/test_sparseLU.cpp
+++ b/bench/spbench/test_sparseLU.cpp
@@ -13,13 +13,14 @@ using namespace Eigen;
 
 int main(int argc, char **args)
 {
-  SparseMatrix<double, ColMajor> A; 
-  typedef SparseMatrix<double, ColMajor>::Index Index;
-  typedef Matrix<double, Dynamic, Dynamic> DenseMatrix;
-  typedef Matrix<double, Dynamic, 1> DenseRhs;
-  VectorXd b, x, tmp;
-//   SparseLU<SparseMatrix<double, ColMajor>, AMDOrdering<int> >   solver;
-  SparseLU<SparseMatrix<double, ColMajor>, COLAMDOrdering<int> >   solver;
+  typedef complex<double> scalar; 
+  SparseMatrix<scalar, ColMajor> A; 
+  typedef SparseMatrix<scalar, ColMajor>::Index Index;
+  typedef Matrix<scalar, Dynamic, Dynamic> DenseMatrix;
+  typedef Matrix<scalar, Dynamic, 1> DenseRhs;
+  Matrix<scalar, Dynamic, 1> b, x, tmp;
+//   SparseLU<SparseMatrix<scalar, ColMajor>, AMDOrdering<int> >   solver;
+  SparseLU<SparseMatrix<scalar, ColMajor>, COLAMDOrdering<int> >   solver;
   ifstream matrix_file; 
   string line;
   int  n;
@@ -36,7 +37,7 @@ int main(int argc, char **args)
   if (iscomplex) { cout<< " Not for complex matrices \n"; return -1; }
   if (isvector) { cout << "The provided file is not a matrix file\n"; return -1;}
   if (sym != 0) { // symmetric matrices, only the lower part is stored
-    SparseMatrix<double, ColMajor> temp; 
+    SparseMatrix<scalar, ColMajor> temp; 
     temp = A;
     A = temp.selfadjointView<Lower>();
   }
@@ -72,8 +73,8 @@ int main(int argc, char **args)
   timer.stop();
   cout << "solve time " << timer.value() << std::endl; 
   /* Check the accuracy */
-  VectorXd tmp2 = b - A*x;
-  double tempNorm = tmp2.norm()/b.norm();
+  Matrix<scalar, Dynamic, 1> tmp2 = b - A*x;
+  scalar tempNorm = tmp2.norm()/b.norm();
   cout << "Relative norm of the computed solution : " << tempNorm <<"\n";
   cout << "Number of nonzeros in the factor : " << solver.nnzL() + solver.nnzU() << std::endl;