Memory management

1 files changed, 93 insertions, 87 deletions

diff --git a/Eigen/src/SparseLU/SparseLU_Memory.h b/Eigen/src/SparseLU/SparseLU_Memory.h
index 730557b63..a92c3bcc4 100644
--- a/Eigen/src/SparseLU/SparseLU_Memory.h
+++ b/Eigen/src/SparseLU/SparseLU_Memory.h
@@ -53,7 +53,7 @@
 namespace internal {
   
 /**
- * \brief  Allocate various working space needed in the numerical factorization phase.
+ * \brief  Allocate various working space failed in the numerical factorization phase.
  * \param m number of rows of the input matrix 
  * \param n number of columns 
  * \param annz number of initial nonzeros in the matrix 
@@ -61,22 +61,21 @@ namespace internal {
  * \param iwork Integer working space 
  * \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; 
- *  FIXME should also return the size of actually allocated when memory allocation failed 
- * NOTE Unlike SuperLU, this routine does not allow the user to provide the size to allocate 
+ * \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 allow the user to provide its own user space
  */
 template <typename ScalarVector,typename IndexVector>
-int SparseLU::LUMemInit(int m, int n, int annz, Scalar *work, Index *iwork, int lwork, int fillratio, GlobalLU_t& Glu)
+int SparseLU::LUMemInit(int m, int n, int annz, ScalarVector& work, IndexVector& iwork, int lwork, int fillratio, GlobalLU_t& Glu)
 {
   typedef typename ScalarVector::Scalar; 
   typedef typename IndexVector::Index; 
   
-  Glu.num_expansions = 0; //No memory expansions so far
-  if (!Glu.expanders)
-    Glu.expanders = new ExpHeader(LU_NBR_MEMTYPE); 
-  
+  int& num_expansions = Glu.num_expansions; //No memory expansions so far
+  num_expansions = 0; 
   // Guess the size for L\U factors 
-  int nzlmax, nzumax, nzlumax; 
+  Index& nzlmax = Glu.nzlmax; 
+  Index& nzumax = Glu.nzumax; 
+  Index& nzlumax = Glu.nzlumax;
   nzumax = nzlumax = m_fillratio * annz; // estimated number of nonzeros in U 
   nzlmax  = std::max(1, m_fill_ratio/4.) * annz; // estimated  nnz in L factor
 
@@ -90,138 +89,145 @@ int SparseLU::LUMemInit(int m, int n, int annz, Scalar *work, Index *iwork, int
   }
   
   // Setup the required space 
-  // NOTE: In SuperLU, there is an option to let the user provide its own space, unlike here.
-  
-  // Allocate Integer pointers for L\U factors
-  Glu.supno = new IndexVector; 
-  Glu.supno->resize(n+1);
-  
-  Glu.xlsub = new IndexVector; 
-  Glu.xlsub->resize(n+1);
   
-  Glu.xlusup = new IndexVector; 
-  Glu.xlusup->resize(n+1);
-  
-  Glu.xusub = new IndexVector; 
-  Glu.xusub->resize(n+1);
+  // First allocate Integer pointers for L\U factors
+  Glu.supno.resize(n+1);
+  Glu.xlsub.resize(n+1);
+  Glu.xlusup.resize(n+1);
+  Glu.xusub.resize(n+1);
 
   // Reserve memory for L/U factors
-  Glu.lusup = new ScalarVector; 
-  Glu.ucol = new ScalarVector; 
-  Glu.lsub = new IndexVector;
-  Glu.usub = new IndexVector; 
-  
-  expand<ScalarVector>(Glu.lusup,nzlumax, LUSUP, 0, 0, Glu); 
-  expand<ScalarVector>(Glu.ucol,nzumax, UCOL, 0, 0, Glu); 
-  expand<IndexVector>(Glu.lsub,nzlmax, LSUB, 0, 0, Glu); 
-  expand<IndexVector>(Glu.usub,nzumax, USUB, 0, 1, Glu); 
+  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, 
   // Should be a try... catch section here 
   while ( !Glu.lusup.size() || !Glu.ucol.size() || !Glu.lsub.size() || !Glu.usub.size())
   {
     //otherwise reduce the estimated size and retry
-//     delete [] Glu.lusup; 
-//     delete [] Glu.ucol;
-//     delete [] Glu.lsub;
-//     delete [] Glu.usub;
-//     
     nzlumax /= 2;
     nzumax /= 2;
     nzlmax /= 2;
-    //FIXME Should be an excpetion here
-    eigen_assert (nzlumax > annz && "Not enough memory to perform factorization");
+    //FIXME Should be an exception here
+    if (nzlumax < annz ) return nzlumax; 
     
-    expand<ScalarVector>(Glu.lsup, nzlumax, LUSUP, 0, 0, Glu); 
-    expand<ScalarVector>(Glu.ucol, nzumax, UCOL, 0, 0, Glu); 
-    expand<IndexVector>(Glu.lsub, nzlmax, LSUB, 0, 0, Glu); 
-    expand<IndexVector>(Glu.usub, nzumax, USUB, 0, 1, Glu); 
+    expand<ScalarVector>(Glu.lsup, nzlumax, 0, 0, Glu); 
+    expand<ScalarVector>(Glu.ucol, nzumax, 0, 0, Glu); 
+    expand<IndexVector>(Glu.lsub, nzlmax, 0, 0, Glu); 
+    expand<IndexVector>(Glu.usub, nzumax, 0, 1, Glu); 
   }
   
   // LUWorkInit : Now, allocate known working storage
   int isize = (2 * m_panel_size + 3 + LU_NO_MARKER) * m + n;
   int dsize = m * m_panel_size + LU_NUM_TEMPV(m, m_panel_size, m_maxsuper, m_rowblk); 
-  iwork = new Index(isize);
-  eigen_assert( (m_iwork != 0) && "Malloc fails for iwork");
-  work = new Scalar(dsize);
-  eigen_assert( (m_work != 0) && "Malloc fails for dwork");
+  iwork.resize(isize); 
+  work.resize(isize); 
   
-  ++Glu.num_expansions;
+  ++num_expansions;
   return 0;
+  
 } // end LuMemInit
 
 /** 
   * Expand the existing storage to accomodate more fill-ins
-  * \param vec Valid pointer to a vector to allocate or expand
-  * \param [in,out]prev_len At input, length from previous call. At output, length of the newly allocated vector
-  * \param type Which part of the memory to expand
-  * \param len_to_copy  Size of the memory to be copied to new store
-  * \param keep_prev  true: use prev_len; Do not expand this vector; false: compute new_len and expand
+  * \param vec Valid pointer to the vector to allocate or expand
+  * \param [in,out]length  At input, contain the current length of the vector that is to be increased. At output, length of the newly allocated vector
+  * \param [in]len_to_copy Current number of elements in the factors
+  * \param keep_prev  true: use length  and do not expand the vector; false: compute new_len and expand
+  * \param [in,out]num_expansions Number of times the memory has been expanded
   */
 template <typename VectorType >
-int  SparseLU::expand(VectorType& vec, int& prev_len, MemType type, int len_to_copy, bool keep_prev, GlobalLU_t& Glu) 
+int  SparseLU::expand(VectorType& vec, int& length, int len_to_copy, bool keep_prev, int& num_expansions) 
 {
   
   float alpha = 1.5; // Ratio of the memory increase 
   int new_len; // New size of the allocated memory
   
-  if(Glu.num_expansions == 0 || keep_prev) 
-    new_len = prev_len; // First time allocate requested
+  if(num_expansions == 0 || keep_prev) 
+    new_len = length ; // First time allocate requested
   else 
-    new_len = alpha * prev_len;
-  
-  // Allocate new space
-//   vec = new VectorType(new_len); 
-  VectorType old_vec(vec); 
-  if ( Glu.num_expansions != 0 ) // The memory has been expanded before
+    new_len = alpha * length ;
+  
+  VectorType old_vec; // Temporary vector to hold the previous values   
+  if (len_to_copy > 0 )
+    old_vec = vec; // old_vec should be of size len_to_copy... to be checked
+  
+  //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
   {
     int tries = 0; 
-    vec.resize(new_len); //expand the current vector 
     if (keep_prev) 
     {
-      if (!vec.size()) return -1 ; // FIXME could throw an exception somehow 
+      if (!vec.size()) return new_len ;
     }
     else 
     {
       while (!vec.size())
       {
-        // Reduce the size and allocate again
-        if ( ++tries > 10) return -1 
+       // Reduce the size and allocate again
+        if ( ++tries > 10) return new_len;  
         alpha = LU_Reduce(alpha); 
-        new_len = alpha * prev_len; 
-        vec->resize(new_len); 
+        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 
-    for (int i = 0; i < old_vec.size(); i++)
-      vec(i) = old_vec(i); 
+    if (len_to_copy > 0)
+      vec.segment(0, len_to_copy) = old_vec;   
   } // end expansion 
-//   expanders[type].mem = vec; 
-//   expanders[type].size = new_len;
-  prev_len = new_len;
-  if(Glu.num_expansions) ++Glu.num_expansions;
+  length  = new_len;
+  if(num_expansions) ++num_expansions;
   return 0; 
 }
 
 /** 
  * \brief Expand the existing storage 
- * 
- * NOTE: The calling sequence of this function is different from that of SuperLU
- * 
- * \return a pointer to the newly allocated space
+ * \param vec vector to expand 
+ * \param [in,out]maxlen On input, previous size of vec (Number of elements to copy ). on output, new size
+ * \param next current number of elements in the vector.
+ * \param Glu Global data structure 
+ * \return 0 on success, > 0 size of the memory allocated so far
  */
-template <typename VectorType>
-VectorType* SparseLU::LUMemXpand(int jcol, int next, MemType mem_type, int& maxlen)
+template <typename IndexVector>
+int SparseLU::LUMemXpand(VectorType& vec, int& maxlen, int next, LU_MemType memtype, LU_GlobalLu_t& Glu)
 {
-  VectorType *newmem; 
+  int failed_size; 
+  int& num_expansions = Glu.num_expansions; 
   if (memtype == USUB)
-    vec = expand<VectorType>(vec, maxlen, mem_type, next, 1);
+     failed_size = expand<IndexVector>(vec, maxlen, next, 1, num_expansions);
   else
-    vec = expand<VectorType>(vec, maxlen, mem_type, next, 0);
-  // FIXME Should be an exception instead of an assert
-  eigen_assert(new_mem.size() && "Can't expand memory"); 
-  
-  return new_mem;
+    failed_size = expand<IndexVector>(vec, maxlen, next, 0, num_expansions);
+
+  if (failed_size)
+    return faileld_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 ; 
   
 }