GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-24 22:04:28 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add exception handler to memory allocation

Browse Source
This commit is contained in:
Desire NUENTSA 2012-07-19 18:03:44 +02:00
parent b0cba2d988
commit 59642da88b
11 changed files with 104 additions and 99 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
Eigen/src/OrderingMethods/Eigen_Colamd.h
View File

@ -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

6
Eigen/src/OrderingMethods/Ordering.h
View File

@ -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:
};

3
Eigen/src/SparseLU/SparseLU.h
View File

@ -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

2
Eigen/src/SparseLU/SparseLU_Coletree.h
View File

@ -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;

1
Eigen/src/SparseLU/SparseLU_Matrix.h
View File

@ -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

157
Eigen/src/SparseLU/SparseLU_Memory.h
View File

@ -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);
//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;
//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;
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 @@ 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;
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 @@ 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 ;
}

1
Eigen/src/SparseLU/SparseLU_column_bmod.h
View File

@ -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);

1
Eigen/src/SparseLU/SparseLU_panel_dfs.h
View File

@ -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;

1
Eigen/src/SparseLU/SparseLU_snode_bmod.h
View File

@ -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;
}

2
bench/spbench/CMakeLists.txt
View File

@ -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})
target_link_libraries (test_sparseLU ${SPARSE_LIBS})

21
bench/spbench/test_sparseLU.cpp
View File

@ -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;