diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index 593ec7e25..996dbf078 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -412,15 +412,14 @@ void SparseLU::factorize(const MatrixType& matrix) // Apply permutation to the L subscripts LU_fixupL(min_mn, m_perm_r, m_Glu); - // Free work space and compress storage iwork and work - // ?? Should it be done automatically by C++ + // Free work space iwork and work //... // Create supernode matrix L m_Lstore.setInfos(m, min_mn, nnzL, Glu.lusup, Glu.xlusup, Glu.lsub, Glu.xlsub, Glu.supno; Glu.xsup); // Create the column major upper sparse matrix U - // Could be great to have the SparseMatrix constructor accepting the CSC matrix pointers - // The Map class can do the job somehow + // ?? Use the MappedSparseMatrix class ?? + new (&m_Ustore) Map > ( m, min_mn, nnzU, Glu.xusub.data(), Glu.usub.data(), Glu.ucol.data() ); m_info = Success; m_factorizationIsOk = ok; } diff --git a/Eigen/src/SparseLU/SparseLU_Matrix.h b/Eigen/src/SparseLU/SparseLU_Matrix.h index 1fe991d1c..01f8784da 100644 --- a/Eigen/src/SparseLU/SparseLU_Matrix.h +++ b/Eigen/src/SparseLU/SparseLU_Matrix.h @@ -53,17 +53,21 @@ class SuperNodalMatrix } SuperNodalMatrix(Index m, Index n, Index nnz, Scalar *nzval, Index* nzval_colptr, Index* rowind, - Index* rowind_colptr, Index* col_to_sup, Index* sup_to_col ):m_row(m),m_col(n),m_nnz(nnz), - m_nzval(nzval),m_nzval_colptr(nzval_colptr),m_rowind(rowind), - m_rowind_colptr(rowind_colptr),m_col_to_sup(col_to_sup),m_sup_to_col(sup_to_col) + Index* rowind_colptr, Index* col_to_sup, Index* sup_to_col ) { - + setInfos(m, n, nnz, nzval, nzval_colptr, rowind, rowind_colptr, col_to_sup, sup_to_col); } ~SuperNodalMatrix() { } + /** + * Set appropriate pointers for the lower triangular supernodal matrix + * These infos are available at the end of the numerical factorization + * FIXME This class will be modified such that it can be use in the course + * of the factorization. + */ void setInfos(Index m, Index n, Index nnz, Scalar *nzval, Index* nzval_colptr, Index* rowind, Index* rowind_colptr, Index* col_to_sup, Index* sup_to_col ) { @@ -78,21 +82,80 @@ class SuperNodalMatrix m_sup_to_col = sup_to_col; } - SuperNodalMatrix(SparseMatrix& mat); - class InnerIterator + /** + * Number of rows + */ + int rows() { - public: - - protected: - - }: + return m_row; + } + + /** + * Number of columns + */ + int cols() + { + return m_col; + } + + /** + * Return the array of nonzero values packed by column + * + * The size is nnz + */ + Scalar* valuePtr() + { + return m_nzval; + } + + /** + * Return the pointers to the beginning of each column in \ref outerIndexPtr() + */ + Index* colIndexPtr() + { + return m_nzval_colptr; + } + + /** + * Return the array of compressed row indices of all supernodes + */ + Index* rowIndex() + { + return m_rowind; + } + /** + * Return the location in \em rowvaluePtr() which starts each column + */ + Index* rowIndexPtr() + { + return m_rowind_colptr; + } + /** + * Return the array of column-to-supernode mapping + */ + Index colToSup() + { + return m_col_to_sup; + } + /** + * Return the array of supernode-to-column mapping + */ + Index supToCol() + { + return m_sup_to_col; + } + + + class InnerIterator; + class SuperNodeIterator; + protected: Index m_row; // Number of rows Index m_col; // Number of columns Index m_nnz; // Number of nonzero values - Index m_nsupper; // Index of the last supernode - Scalar* m_nzval; //array of nonzero values packed by (supernode ??) column + 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 Index* m_rowind; // Array of compressed row indices of rectangular supernodes Index* m_rowind_colptr; //rowind_colptr[j] stores the location in rowind[] which starts column j @@ -100,11 +163,90 @@ class SuperNodalMatrix Index *m_sup_to_col; //sup_to_col[s] points to the starting column of the s-th supernode private : - SuperNodalMatrix(SparseMatrix& ) {} }; -SuperNodalMatrix::SuperNodalMatrix(SparseMatrix& mat) +/** + * \brief InnerIterator class to iterate over nonzero values in the triangular supernodal matrix + * + */ +template +class SuperNodalMatrix::InnerIterator { + public: + InnerIterator(const SuperNodalMatrix& mat, Index outer) + : m_matrix(mat), + m_outer(outer), + m_idval(mat.colIndexPtr()[outer]), + m_startval(m_idval), + m_endval(mat.colIndexPtr()[outer+1]) + m_idrow(mat.rowIndexPtr()[outer]), + m_startidrow(m_idrow), + m_endidrow(mat.rowIndexPtr()[outer+1]) + {} + inline InnerIterator& operator++() + { + m_idval++; + m_idrow++ ; + return *this; + } + inline Scalar value() const { return m_matrix.valuePtr()[m_idval]; } + + inline Scalar& valueRef() { return const_cast(m_matrix.valuePtr()[m_idval]; } + + inline Index index() const { return m_matrix.rowIndex()[m_idrow]; } + inline Index row() const { return index(); } + inline Index col() const { return m_outer; } + + inline Index supIndex() const { return m_matrix.colToSup()[m_outer]; } + + inline operator bool() const + { + return ( (m_idval < m_endval) && (m_idval > m_startval) && + (m_idrow < m_endidrow) && (m_idrow > m_startidrow) ); + } + + protected: + const SuperNodalMatrix& m_matrix; // Supernodal lower triangular matrix + const Index m_outer; // Current column + Index m_idval; //Index to browse the values in the current column + const Index m_startval; // Start of the column value + const Index m_endval; // End of the column value + Index m_idrow; //Index to browse the row indices + const Index m_startidrow; // Start of the row indices of the current column value + const Index m_endidrow; // End of the row indices of the current column value +}; +/** + * \brief Iterator class to iterate over nonzeros Supernodes in the triangular supernodal matrix + * + * The final goal is to use this class when dealing with supernodes during numerical factorization + */ +template +class SuperNodalMatrix::SuperNodeIterator +{ + public: + SuperNodeIterator(const SuperNodalMatrix& mat) + { + + } + SuperNodeIterator(const SuperNodalMatrix& mat, Index supno) + { + + } + + /* + * Available Methods : + * Browse all supernodes (operator ++ ) + * Number of supernodes + * Columns of the current supernode + * triangular matrix of the current supernode + * rectangular part of the current supernode + */ + protected: + const SuperNodalMatrix& m_matrix; // Supernodal lower triangular matrix + Index m_idsup; // Index to browse all supernodes + const Index m_nsuper; // Number of all supernodes + Index m_startidsup; + Index m_endidsup; -} +}; #endif \ No newline at end of file diff --git a/Eigen/src/SparseLU/SparseLU_Memory.h b/Eigen/src/SparseLU/SparseLU_Memory.h index 91b24fa67..730557b63 100644 --- a/Eigen/src/SparseLU/SparseLU_Memory.h +++ b/Eigen/src/SparseLU/SparseLU_Memory.h @@ -52,169 +52,155 @@ + (w + 1) * m * sizeof(Scalar) namespace internal { -/* Allocate various working space needed in the numerical factorization phase. - * m_work : space fot the output data structures (lwork is the size) - * m_Glu: persistent data to facilitate multiple factors : is it really necessary ?? +/** + * \brief Allocate various working space needed 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 + * \param work scalar working space needed by all factor routines + * \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 - * nor it return an estimated amount of space required. - * - * Useful variables : - * - m_fillratio : Ratio of fill expected - * - lwork = -1 : return an estimated size of the required memory - * = 0 : Estimate and allocate the memory */ -template -int SparseLU::LUMemInit(int lwork) +template +int SparseLU::LUMemInit(int m, int n, int annz, Scalar *work, Index *iwork, int lwork, int fillratio, GlobalLU_t& Glu) { - int iword = sizeof(Index); - int dword = sizeof(Scalar); - int n = m_Glu.n = m_mat.cols(); - int m = m_mat.rows(); - m_Glu.num_expansions = 0; // No memory expansions so far ?? + 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); + + // Guess the size for L\U factors + int nzlmax, nzumax, 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 + + // Return the estimated size to the user if necessary int estimated_size; - - - if (!m_Glu.expanders) - m_Glu.expanders = new ExpHeader(NO_MEMTYPE); - - if (m_fact_t != SamePattern_SameRowPerm) // Create space for a new factorization - { - // Guess the size for L\U factors - int annz = m_mat.nonZeros(); - int nzlmax, nzumax, nzlumax; - nzumax = nzlumax = m_fillratio * annz; // ??? - nzlmax = std::max(1, m_fill_ratio/4.) * annz; //??? - - // Return the estimated size to the user if necessary - if (lwork == IND_EMPTY) - { - estimated_size = LU_GluIntArray(n) * iword + LU_TempSpace(m, m_panel_size) - + (nzlmax + nzumax) * iword + (nzlumax+nzumax) * dword + n); - return estimated_size; - } - - // Setup the required space - // NOTE: In SuperLU, there is an option to let the user provide its own space. - - // Allocate Integer pointers for L\U factors.resize(n+1); - m_Glu.supno.resize(n+1); - m_Glu.xlsub.resize(n+1); - m_Glu.xlusup.resize(n+1); - m_Glu.xusub.resize(n+1); - - // Reserve memory for L/U factors - m_Glu.lusup = internal::expand(nzlumax, LUSUP, 0, 0, m_Glu); - m_Glu.ucol = internal::expand(nzumax, UCOL, 0, 0, m_Glu); - m_Glu.lsub = internal::expand(nzlmax, LSUB, 0, 0, m_Glu); - m_Glu.usub = internal::expand(nzumax, USUB, 0, 1, m_Glu); - - // Check if the memory is correctly allocated, - while ( !m_Glu.lusup || !m_Glu.ucol || !m_Glu.lsub || !m_Glu.usub) - { - //otherwise reduce the estimated size and retry - delete [] m_Glu.lusup; - delete [] m_Glu.ucol; - delete [] m_Glu.lsub; - delete [] m_Glu.usub; - - nzlumax /= 2; - nzumax /= 2; - nzlmax /= 2; - eigen_assert (nzlumax > annz && "Not enough memory to perform factorization"); - - m_Glu.lusup = internal::expand(nzlumax, LUSUP, 0, 0, m_Glu); - m_Glu.ucol = internal::expand(nzumax, UCOL, 0, 0, m_Glu); - m_Glu.lsub = internal::expand(nzlmax, LSUB, 0, 0, m_Glu); - m_Glu.usub = internal::expand(nzumax, USUB, 0, 1, m_Glu); - } - } - else // m_fact == SamePattern_SameRowPerm; + if (lwork == IND_EMPTY) { - if (lwork == IND_EMPTY) - { - estimated_size = LU_GluIntArray(n) * iword + LU_TempSpace(m, m_panel_size) - + (Glu.nzlmax + Glu.nzumax) * iword + (Glu.nzlumax+Glu.nzumax) * dword + n); - return estimated_size; - } - // Use existing space from previous factorization - // Unlike in SuperLU, this should not be necessary here since m_Glu is persistent as a member of the class - m_Glu.xsup = m_Lstore.sup_to_col; - m_Glu.supno = m_Lstore.col_to_sup; - m_Glu.xlsub = m_Lstore.rowind_colptr; - m_Glu.xlusup = m_Lstore.nzval_colptr; - xusub = m_Ustore.outerIndexPtr(); + estimated_size = LU_GluIntArray(n) * sizeof(Index) + LU_TempSpace(m, m_panel_size) + + (nzlmax + nzumax) * sizeof(Index) + (nzlumax+nzumax) * sizeof(Scalar) + n); + return estimated_size; + } + + // 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); + + // Reserve memory for L/U factors + Glu.lusup = new ScalarVector; + Glu.ucol = new ScalarVector; + Glu.lsub = new IndexVector; + Glu.usub = new IndexVector; + + expand(Glu.lusup,nzlumax, LUSUP, 0, 0, Glu); + expand(Glu.ucol,nzumax, UCOL, 0, 0, Glu); + expand(Glu.lsub,nzlmax, LSUB, 0, 0, Glu); + expand(Glu.usub,nzumax, USUB, 0, 1, Glu); + + // 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"); - m_Glu.expanders[LSUB].size = m_Glu.nzlmax; // Maximum value from previous factorization - m_Glu.expanders[LUSUP].size = m_Glu.nzlumax; - m_Glu.expanders[USUB].size = GLu.nzumax; - m_Glu.expanders[UCOL].size = m_Glu.nzumax; - m_Glu.lsub = GLu.expanders[LSUB].mem = m_Lstore.rowind; - m_Glu.lusup = GLu.expanders[LUSUP].mem = m_Lstore.nzval; - GLu.usub = m_Glu.expanders[USUB].mem = m_Ustore.InnerIndexPtr(); - m_Glu.ucol = m_Glu.expanders[UCOL].mem = m_Ustore.valuePtr(); + expand(Glu.lsup, nzlumax, LUSUP, 0, 0, Glu); + expand(Glu.ucol, nzumax, UCOL, 0, 0, Glu); + expand(Glu.lsub, nzlmax, LSUB, 0, 0, Glu); + expand(Glu.usub, nzumax, USUB, 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); - m_iwork = new Index(isize); + iwork = new Index(isize); eigen_assert( (m_iwork != 0) && "Malloc fails for iwork"); - m_work = new Scalar(dsize); + work = new Scalar(dsize); eigen_assert( (m_work != 0) && "Malloc fails for dwork"); - ++m_Glu.num_expansions; + ++Glu.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 */ -template -DestType* SparseLU::expand(int& prev_len, // Length from previous call - MemType type, // Which part of the memory to expand - int len_to_copy, // Size of the memory to be copied to new store - int keep_prev) // = 1: use prev_len; Do not expand this vector - // = 0: compute new_len to expand) +template +int SparseLU::expand(VectorType& vec, int& prev_len, MemType type, int len_to_copy, bool keep_prev, GlobalLU_t& Glu) { float alpha = 1.5; // Ratio of the memory increase int new_len; // New size of the allocated memory - if(m_Glu.num_expansions == 0 || keep_prev) - new_len = prev_len; + + if(Glu.num_expansions == 0 || keep_prev) + new_len = prev_len; // First time allocate requested else new_len = alpha * prev_len; - // Allocate new space - DestType *new_mem, *old_mem; - new_mem = new DestType(new_len); - if ( m_Glu.num_expansions != 0 ) // The memory has been expanded before + // Allocate new space +// vec = new VectorType(new_len); + VectorType old_vec(vec); + if ( Glu.num_expansions != 0 ) // The memory has been expanded before { int tries = 0; + vec.resize(new_len); //expand the current vector if (keep_prev) { - if (!new_mem) return 0; + if (!vec.size()) return -1 ; // FIXME could throw an exception somehow } else { - while ( !new_mem) + while (!vec.size()) { // Reduce the size and allocate again - if ( ++tries > 10) return 0; + if ( ++tries > 10) return -1 alpha = LU_Reduce(alpha); new_len = alpha * prev_len; - new_mem = new DestType(new_len); + vec->resize(new_len); } - } // keep_prev + } // end allocation //Copy the previous values to the newly allocated space - ExpHeader* expanders = m_Glu.expanders; - std::memcpy(new_mem, expanders[type].mem, len_to_copy); - delete [] expanders[type].mem; - } - expanders[type].mem = new_mem; - expanders[type].size = new_len; + for (int i = 0; i < old_vec.size(); i++) + vec(i) = old_vec(i); + } // end expansion +// expanders[type].mem = vec; +// expanders[type].size = new_len; prev_len = new_len; - if(m_Glu.num_expansions) ++m_Glu.num_expansions; - return expanders[type].mem; + if(Glu.num_expansions) ++Glu.num_expansions; + return 0; } /** @@ -224,15 +210,16 @@ DestType* SparseLU::expand(int& prev_len, // Length from previous call * * \return a pointer to the newly allocated space */ -template -DestType* SparseLU::LUMemXpand(int jcol, int next, MemType mem_type, int& maxlen) +template +VectorType* SparseLU::LUMemXpand(int jcol, int next, MemType mem_type, int& maxlen) { - DestType *newmem; + VectorType *newmem; if (memtype == USUB) - new_mem = expand(maxlen, mem_type, next, 1); + vec = expand(vec, maxlen, mem_type, next, 1); else - new_mem = expand(maxlen, mem_type, next, 0); - eigen_assert(new_mem && "Can't expand memory"); // FIXME Should be an exception + vec = expand(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; diff --git a/Eigen/src/SparseLU/SparseLU_Structs.h b/Eigen/src/SparseLU/SparseLU_Structs.h index 72e1db343..e680eaa21 100644 --- a/Eigen/src/SparseLU/SparseLU_Structs.h +++ b/Eigen/src/SparseLU/SparseLU_Structs.h @@ -23,7 +23,7 @@ // Eigen. If not, see . /* - * NOTE: Part of this file is the modified version of files slu_[s,d,c,z]defs.h + * NOTE: This file comes from a partly modified version of files slu_[s,d,c,z]defs.h * -- SuperLU routine (version 4.1) -- * Univ. of California Berkeley, Xerox Palo Alto Research Center, * and Lawrence Berkeley National Lab. @@ -84,36 +84,39 @@ #define EIGEN_LU_STRUCTS namespace Eigen { -#define NO_MEMTYPE 4 /* 0: lusup +#define LU_NBR_MEMTYPE 4 /* 0: lusup 1: ucol 2: lsub 3: usub */ -typedef enum {NATURAL, MMD_ATA, MMD_AT_PLUS_A, COLAMD, MY_PREMC} colperm_t; -typedef enum {DOFACT, SamePattern, SamePattern_SameRowPerm, Factored} fact_t; +typedef enum {NATURAL, MMD_ATA, MMD_AT_PLUS_A, COLAMD, MY_PERMC} colperm_t; +typedef enum {DOFACT, SamePattern, Factored} fact_t; typedef enum {LUSUP, UCOL, LSUB, USUB, LLVL, ULVL} MemType; /** Headers for dynamically managed memory - \tparam BaseType can be int, real scalar or complex scalar*/ -template + \tparam IndexVectorType can be int, real scalar or complex scalar*/ +template struct ExpHeader { int size; // Length of the memory that has been used */ - BaseType *mem; + VectorType *mem; // Save the current pointer of the newly allocated memory } ExpHeader; -template +template struct { - VectorXi xsup; // supernode and column mapping - VectorXi supno; // Supernode number corresponding to this column - VectorXi lsub; // Compressed L subscripts of rectangular supernodes - VectorXi xlsub; // xlsub(j) points to the starting location of the j-th column in lsub - VectorXi xlusup; - VectorXi xusub; - VectorType lusup; // L supernodes - VectorType ucol; // U columns + IndexVector* xsup; //First supernode column ... xsup(s) points to the beginning of the s-th supernode + IndexVector* supno; // Supernode number corresponding to this column (column to supernode mapping) + ScalarVector* lusup; // nonzero values of L ordered by columns + IndexVector* lsub; // Compressed row indices of L rectangular supernodes. + IndexVector* xlusup; // pointers to the beginning of each column in lusup + IndexVector* xlsub; // pointers to the beginning of each column in lsub Index nzlmax; // Current max size of lsub - Index nzumax; // Current max size of ucol Index nzlumax; // Current max size of lusup + + ScalarVector* ucol; // nonzero values of U ordered by columns + IndexVector* usub; // row indices of U columns in ucol + IndexVector* xusub; // Pointers to the beginning of each column of U in ucol + Index nzumax; // Current max size of ucol Index n; // Number of columns in the matrix + int num_expansions; ExpHeader *expanders; // Array of pointers to 4 types of memory } GlobalLU_t;