Add exception handler to memory allocation
diff --git a/Eigen/src/OrderingMethods/Eigen_Colamd.h b/Eigen/src/OrderingMethods/Eigen_Colamd.h index 39701d0..0af137d 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 @@ 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 @@ } #endif /* NDEBUG */ - #endif
diff --git a/Eigen/src/OrderingMethods/Ordering.h b/Eigen/src/OrderingMethods/Ordering.h index cbd2e5d..47cd6f1 100644 --- a/Eigen/src/OrderingMethods/Ordering.h +++ b/Eigen/src/OrderingMethods/Ordering.h
@@ -124,9 +124,6 @@ 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 @@ 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 bb1decc..25fad0f 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h
@@ -339,9 +339,6 @@ //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 1329d38..142f499 100644 --- a/Eigen/src/SparseLU/SparseLU_Coletree.h +++ b/Eigen/src/SparseLU/SparseLU_Coletree.h
@@ -94,7 +94,6 @@ 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 @@ 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 90a0f27..9f2dcaa 100644 --- a/Eigen/src/SparseLU/SparseLU_Matrix.h +++ b/Eigen/src/SparseLU/SparseLU_Matrix.h
@@ -192,7 +192,6 @@ 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 a170791..7a2ab93 100644 --- a/Eigen/src/SparseLU/SparseLU_Memory.h +++ b/Eigen/src/SparseLU/SparseLU_Memory.h
@@ -78,41 +78,82 @@ 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); + //Allocate or expand the current vector + try + { + vec.resize(new_len); + } + catch(std::bad_alloc& ) + { + if ( !num_expansions ) + { + // 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 + { + // Reduce the size and increase again + int tries = 0; // Number of attempts + do + { + alpha = LU_Reduce(alpha); + new_len = alpha * length ; + 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 - length = new_len; - if(num_expansions) ++num_expansions; - return 0; +// 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 @@ * \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 @@ 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; + 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; + } - 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 @@ 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 3042eb5..0078772 100644 --- a/Eigen/src/SparseLU/SparseLU_column_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_column_bmod.h
@@ -133,7 +133,6 @@ // 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 908ee67..70ea0f5 100644 --- a/Eigen/src/SparseLU/SparseLU_panel_dfs.h +++ b/Eigen/src/SparseLU/SparseLU_panel_dfs.h
@@ -92,7 +92,6 @@ 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 44438d0..fc8042f 100644 --- a/Eigen/src/SparseLU/SparseLU_snode_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_snode_bmod.h
@@ -80,7 +80,6 @@ // 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 4b3c6f8..a093cc5 100644 --- a/bench/spbench/CMakeLists.txt +++ b/bench/spbench/CMakeLists.txt
@@ -67,4 +67,4 @@ 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 6fbf034..31273ad 100644 --- a/bench/spbench/test_sparseLU.cpp +++ b/bench/spbench/test_sparseLU.cpp
@@ -13,13 +13,14 @@ 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 @@ 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 @@ 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;