diff options
-rw-r--r-- | Eigen/src/OrderingMethods/Eigen_Colamd.h | 8 | ||||
-rw-r--r-- | Eigen/src/OrderingMethods/Ordering.h | 6 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU.h | 3 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_Coletree.h | 2 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_Matrix.h | 1 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_Memory.h | 147 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_column_bmod.h | 1 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_panel_dfs.h | 1 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_snode_bmod.h | 1 | ||||
-rw-r--r-- | bench/spbench/CMakeLists.txt | 2 | ||||
-rw-r--r-- | bench/spbench/test_sparseLU.cpp | 21 |
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; |