From a01371548dc66ee8cbfac8effd5f560bf5d5697a Mon Sep 17 00:00:00 2001 From: Desire NUENTSA Date: Tue, 25 Sep 2012 09:53:40 +0200 Subject: Define sparseLU functions as static --- Eigen/src/SparseLU/SparseLU.h | 97 ++++++++++---------------- Eigen/src/SparseLU/SparseLU_Coletree.h | 16 ++--- Eigen/src/SparseLU/SparseLU_Memory.h | 15 ++-- Eigen/src/SparseLU/SparseLU_Utils.h | 8 +-- Eigen/src/SparseLU/SparseLU_column_bmod.h | 9 +-- Eigen/src/SparseLU/SparseLU_column_dfs.h | 17 +++-- Eigen/src/SparseLU/SparseLU_copy_to_ucol.h | 8 +-- Eigen/src/SparseLU/SparseLU_heap_relax_snode.h | 4 +- Eigen/src/SparseLU/SparseLU_kernel_bmod.h | 2 +- Eigen/src/SparseLU/SparseLU_panel_bmod.h | 12 ++-- Eigen/src/SparseLU/SparseLU_panel_dfs.h | 14 ++-- Eigen/src/SparseLU/SparseLU_pivotL.h | 7 +- Eigen/src/SparseLU/SparseLU_pruneL.h | 7 +- Eigen/src/SparseLU/SparseLU_relax_snode.h | 4 +- Eigen/src/SparseLU/SparseLU_snode_bmod.h | 8 +-- Eigen/src/SparseLU/SparseLU_snode_dfs.h | 81 +++++++++++---------- 16 files changed, 130 insertions(+), 179 deletions(-) (limited to 'Eigen/src/SparseLU') diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index 77df091c3..f5d15ec6b 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -18,6 +18,8 @@ namespace Eigen { #include "SparseLU_Structs.h" #include "SparseLU_Matrix.h" +// Base structure containing all the factorization routines +#include "SparseLUBase.h" /** * \ingroup SparseLU_Module * \brief Sparse supernodal LU factorization for general matrices @@ -40,6 +42,7 @@ class SparseLU typedef Matrix ScalarVector; typedef Matrix IndexVector; typedef PermutationMatrix PermutationType; + public: SparseLU():m_isInitialized(true),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0) { @@ -58,6 +61,7 @@ class SparseLU void analyzePattern (const MatrixType& matrix); void factorize (const MatrixType& matrix); + void simplicialfactorize(const MatrixType& matrix); /** * Compute the symbolic and numeric factorization of the input sparse matrix. @@ -224,8 +228,7 @@ class SparseLU PermutationType m_perm_r ; // Row permutation IndexVector m_etree; // Column elimination tree - LU_GlobalLU_t m_glu; // persistent data to facilitate multiple factors - // FIXME All fields of this struct can be defined separately as class members + LU_GlobalLU_t m_glu; // SuperLU/SparseLU options bool m_symmetricmode; @@ -243,7 +246,6 @@ class SparseLU // Functions needed by the anaysis phase -#include "SparseLU_Coletree.h" /** * Compute the column permutation to minimize the fill-in (file amd.c ) * @@ -262,9 +264,6 @@ void SparseLU::analyzePattern(const MatrixType& mat) OrderingType ord; ord(mat,m_perm_c); - //FIXME Check the right semantic behind m_perm_c - // that is, column j of mat goes to column m_perm_c(j) of mat * m_perm_c; - // Apply the permutation to the column of the input matrix // m_mat = mat * m_perm_c.inverse(); //FIXME It should be less expensive here to permute only the structural pattern of the matrix @@ -282,13 +281,13 @@ void SparseLU::analyzePattern(const MatrixType& mat) // Compute the column elimination tree of the permuted matrix /*if (m_etree.size() == 0) */m_etree.resize(m_mat.cols()); - LU_sp_coletree(m_mat, m_etree); + SparseLUBase::LU_sp_coletree(m_mat, m_etree); // In symmetric mode, do not do postorder here if (!m_symmetricmode) { IndexVector post, iwork; // Post order etree - LU_TreePostorder(m_mat.cols(), m_etree, post); + SparseLUBase::LU_TreePostorder(m_mat.cols(), m_etree, post); // Renumber etree in postorder @@ -310,21 +309,7 @@ void SparseLU::analyzePattern(const MatrixType& mat) m_analysisIsOk = true; } -// Functions needed by the numerical factorization phase -#include "SparseLU_Memory.h" -#include "SparseLU_heap_relax_snode.h" -#include "SparseLU_relax_snode.h" -#include "SparseLU_snode_dfs.h" -#include "SparseLU_snode_bmod.h" -#include "SparseLU_pivotL.h" -#include "SparseLU_panel_dfs.h" -#include "SparseLU_kernel_bmod.h" -#include "SparseLU_panel_bmod.h" -#include "SparseLU_column_dfs.h" -#include "SparseLU_column_bmod.h" -#include "SparseLU_copy_to_ucol.h" -#include "SparseLU_pruneL.h" -#include "SparseLU_Utils.h" +// Functions needed by the numerical factorization phase /** @@ -370,7 +355,7 @@ void SparseLU::factorize(const MatrixType& matrix) int maxpanel = m_perfv.panel_size * m; // Allocate working storage common to the factor routines int lwork = 0; - int info = LUMemInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu); + int info = SparseLUBase::LUMemInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu); if (info) { std::cerr << "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ; @@ -402,25 +387,17 @@ void SparseLU::factorize(const MatrixType& matrix) // Identify initial relaxed snodes IndexVector relax_end(n); if ( m_symmetricmode == true ) - LU_heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end); + SparseLUBase::LU_heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end); else - LU_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end); + SparseLUBase::LU_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end); m_perm_r.resize(m); m_perm_r.indices().setConstant(-1); marker.setConstant(-1); - IndexVector& xsup = m_glu.xsup; - IndexVector& supno = m_glu.supno; - IndexVector& xlsub = m_glu.xlsub; - IndexVector& xlusup = m_glu.xlusup; - IndexVector& xusub = m_glu.xusub; - ScalarVector& lusup = m_glu.lusup; - Index& nzlumax = m_glu.nzlumax; - - supno(0) = IND_EMPTY; xsup.setConstant(0); - xsup(0) = xlsub(0) = xusub(0) = xlusup(0) = Index(0); + m_glu.supno(0) = IND_EMPTY; m_glu.xsup.setConstant(0); + m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0); // Work on one 'panel' at a time. A panel is one of the following : // (a) a relaxed supernode at the bottom of the etree, or @@ -441,7 +418,7 @@ void SparseLU::factorize(const MatrixType& matrix) // Factorize the relaxed supernode(jcol:kcol) // First, determine the union of the row structure of the snode - info = LU_snode_dfs(jcol, kcol, m_mat, xprune, marker, m_glu); + info = SparseLUBase::LU_snode_dfs(jcol, kcol, m_mat, xprune, marker, m_glu); if ( info ) { std::cerr << "MEMORY ALLOCATION FAILED IN SNODE_DFS() \n"; @@ -449,15 +426,15 @@ void SparseLU::factorize(const MatrixType& matrix) m_factorizationIsOk = false; return; } - nextu = xusub(jcol); //starting location of column jcol in ucol - nextlu = xlusup(jcol); //Starting location of column jcol in lusup (rectangular supernodes) - jsupno = supno(jcol); // Supernode number which column jcol belongs to - fsupc = xsup(jsupno); //First column number of the current supernode - new_next = nextlu + (xlsub(fsupc+1)-xlsub(fsupc)) * (kcol - jcol + 1); + nextu = m_glu.xusub(jcol); //starting location of column jcol in ucol + nextlu = m_glu.xlusup(jcol); //Starting location of column jcol in lusup (rectangular supernodes) + jsupno = m_glu.supno(jcol); // Supernode number which column jcol belongs to + fsupc = m_glu.xsup(jsupno); //First column number of the current supernode + new_next = nextlu + (m_glu.xlsub(fsupc+1)-m_glu.xlsub(fsupc)) * (kcol - jcol + 1); int mem; - while (new_next > nzlumax ) + while (new_next > m_glu.nzlumax ) { - mem = LUMemXpand(lusup, nzlumax, nextlu, LUSUP, m_glu.num_expansions); + mem = SparseLUBase::LUMemXpand(m_glu.lusup, m_glu.nzlumax, nextlu, LUSUP, m_glu.num_expansions); if (mem) { std::cerr << "MEMORY ALLOCATION FAILED FOR L FACTOR \n"; @@ -468,16 +445,16 @@ void SparseLU::factorize(const MatrixType& matrix) // Now, left-looking factorize each column within the snode for (icol = jcol; icol<=kcol; icol++){ - xusub(icol+1) = nextu; + m_glu.xusub(icol+1) = nextu; // Scatter into SPA dense(*) for (typename MatrixType::InnerIterator it(m_mat, icol); it; ++it) dense(it.row()) = it.value(); // Numeric update within the snode - LU_snode_bmod(icol, fsupc, dense, m_glu); + SparseLUBase::LU_snode_bmod(icol, fsupc, dense, m_glu); // Eliminate the current column - info = LU_pivotL(icol, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu); + info = SparseLUBase::LU_pivotL(icol, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu); if ( info ) { m_info = NumericalIssue; @@ -505,10 +482,10 @@ void SparseLU::factorize(const MatrixType& matrix) panel_size = n - jcol; // Symbolic outer factorization on a panel of columns - LU_panel_dfs(m, panel_size, jcol, m_mat, m_perm_r.indices(), nseg1, dense, panel_lsub, segrep, repfnz, xprune, marker, parent, xplore, m_glu); + SparseLUBase::LU_panel_dfs(m, panel_size, jcol, m_mat, m_perm_r.indices(), nseg1, dense, panel_lsub, segrep, repfnz, xprune, marker, parent, xplore, m_glu); // Numeric sup-panel updates in topological order - LU_panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_perfv, m_glu); + SparseLUBase::LU_panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_perfv, m_glu); // Sparse LU within the panel, and below the panel diagonal for ( jj = jcol; jj< jcol + panel_size; jj++) @@ -519,7 +496,7 @@ void SparseLU::factorize(const MatrixType& matrix) //Depth-first-search for the current column VectorBlock panel_lsubk(panel_lsub, k, m); VectorBlock repfnz_k(repfnz, k, m); - info = LU_column_dfs(m, jj, m_perm_r.indices(), m_perfv.maxsuper, nseg, panel_lsubk, segrep, repfnz_k, xprune, marker, parent, xplore, m_glu); + info = SparseLUBase::LU_column_dfs(m, jj, m_perm_r.indices(), m_perfv.maxsuper, nseg, panel_lsubk, segrep, repfnz_k, xprune, marker, parent, xplore, m_glu); if ( info ) { std::cerr << "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() \n"; @@ -530,7 +507,7 @@ void SparseLU::factorize(const MatrixType& matrix) // Numeric updates to this column VectorBlock dense_k(dense, k, m); VectorBlock segrep_k(segrep, nseg1, m-nseg1); - info = LU_column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu); + info = SparseLUBase::LU_column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu); if ( info ) { std::cerr << "UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() \n"; @@ -540,7 +517,7 @@ void SparseLU::factorize(const MatrixType& matrix) } // Copy the U-segments to ucol(*) - info = LU_copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu); + info = SparseLUBase::LU_copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu); if ( info ) { std::cerr << "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() \n"; @@ -550,7 +527,7 @@ void SparseLU::factorize(const MatrixType& matrix) } // Form the L-segment - info = LU_pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu); + info = SparseLUBase::LU_pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu); if ( info ) { std::cerr<< "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT " << info <::factorize(const MatrixType& matrix) } // Prune columns (0:jj-1) using column jj - LU_pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu); + SparseLUBase::LU_pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu); // Reset repfnz for this column for (i = 0; i < nseg; i++) @@ -574,11 +551,9 @@ void SparseLU::factorize(const MatrixType& matrix) } // end for -- end elimination // Count the number of nonzeros in factors - LU_countnz(n, m_nnzL, m_nnzU, m_glu); + SparseLUBase::LU_countnz(n, m_nnzL, m_nnzU, m_glu); // Apply permutation to the L subscripts - LU_fixupL(n, m_perm_r.indices(), m_glu); - - + SparseLUBase::LU_fixupL(n, m_perm_r.indices(), m_glu); // Create supernode matrix L m_Lstore.setInfos(m, n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup); @@ -589,7 +564,7 @@ void SparseLU::factorize(const MatrixType& matrix) m_factorizationIsOk = true; } - +// #include "SparseLU_simplicialfactorize.h" namespace internal { template @@ -607,7 +582,5 @@ struct solve_retval, Rhs> } // end namespace internal - - } // End namespace Eigen -#endif \ No newline at end of file +#endif diff --git a/Eigen/src/SparseLU/SparseLU_Coletree.h b/Eigen/src/SparseLU/SparseLU_Coletree.h index 964f5e433..bb4067a45 100644 --- a/Eigen/src/SparseLU/SparseLU_Coletree.h +++ b/Eigen/src/SparseLU/SparseLU_Coletree.h @@ -31,8 +31,8 @@ #ifndef SPARSELU_COLETREE_H #define SPARSELU_COLETREE_H /** Find the root of the tree/set containing the vertex i : Use Path halving */ -template -int etree_find (int i, IndexVector& pp) +template< typename Scalar,typename Index> +int SparseLUBase::etree_find (int i, IndexVector& pp) { int p = pp(i); // Parent int gp = pp(p); // Grand parent @@ -50,8 +50,8 @@ int etree_find (int i, IndexVector& pp) * NOTE : The matrix is supposed to be in column-major format. * */ -template -int LU_sp_coletree(const MatrixType& mat, IndexVector& parent) +template +int SparseLUBase::LU_sp_coletree(const MatrixType& mat, IndexVector& parent) { int nc = mat.cols(); // Number of columns int nr = mat.rows(); // Number of rows @@ -106,8 +106,8 @@ int LU_sp_coletree(const MatrixType& mat, IndexVector& parent) * Depth-first search from vertex n. No recursion. * This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France. */ -template -void LU_nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, int postnum) +template +void SparseLUBase::LU_nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, int postnum) { int current = n, first, next; while (postnum != n) @@ -152,8 +152,8 @@ void LU_nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVecto * \param parent Input tree * \param post postordered tree */ -template -void LU_TreePostorder(int n, IndexVector& parent, IndexVector& post) +template +void SparseLUBase::LU_TreePostorder(int n, IndexVector& parent, IndexVector& post) { IndexVector first_kid, next_kid; // Linked list of children int postnum; diff --git a/Eigen/src/SparseLU/SparseLU_Memory.h b/Eigen/src/SparseLU/SparseLU_Memory.h index 48b36f5b4..0396ab61f 100644 --- a/Eigen/src/SparseLU/SparseLU_Memory.h +++ b/Eigen/src/SparseLU/SparseLU_Memory.h @@ -49,8 +49,9 @@ * \param keep_prev 1: use length and do not expand the vector; 0: compute new_len and expand * \param [in,out]num_expansions Number of times the memory has been expanded */ -template -int expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_expansions) +template +template +int SparseLUBase::expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_expansions) { float alpha = 1.5; // Ratio of the memory increase @@ -122,12 +123,9 @@ int expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_ex * \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 -int LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size, LU_GlobalLU_t& glu) +template +int SparseLUBase::LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size, GlobalLU_t& glu) { - typedef typename ScalarVector::Scalar Scalar; - typedef typename IndexVector::Scalar Index; - int& num_expansions = glu.num_expansions; //No memory expansions so far num_expansions = 0; glu.nzumax = glu.nzlumax = std::max(fillratio * annz, m*n); // estimated number of nonzeros in U @@ -187,8 +185,9 @@ int LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size, * \param glu Global data structure * \return 0 on success, > 0 size of the memory allocated so far */ +template template -int LUMemXpand(VectorType& vec, int& maxlen, int nbElts, LU_MemType memtype, int& num_expansions) +int SparseLUBase::LUMemXpand(VectorType& vec, int& maxlen, int nbElts, LU_MemType memtype, int& num_expansions) { int failed_size; if (memtype == USUB) diff --git a/Eigen/src/SparseLU/SparseLU_Utils.h b/Eigen/src/SparseLU/SparseLU_Utils.h index 316b09ab0..b13930dbb 100644 --- a/Eigen/src/SparseLU/SparseLU_Utils.h +++ b/Eigen/src/SparseLU/SparseLU_Utils.h @@ -15,8 +15,8 @@ /** * \brief Count Nonzero elements in the factors */ -template -void LU_countnz(const int n, int& nnzL, int& nnzU, LU_GlobalLU_t& glu) +template +void SparseLUBase::LU_countnz(const int n, int& nnzL, int& nnzU, GlobalLU_t& glu) { nnzL = 0; nnzU = (glu.xusub)(n); @@ -46,8 +46,8 @@ void LU_countnz(const int n, int& nnzL, int& nnzU, LU_GlobalLU_t -void LU_fixupL(const int n, const IndexVector& perm_r, LU_GlobalLU_t& glu) +template +void SparseLUBase::LU_fixupL(const int n, const IndexVector& perm_r, GlobalLU_t& glu) { int fsupc, i, j, k, jstart; diff --git a/Eigen/src/SparseLU/SparseLU_column_bmod.h b/Eigen/src/SparseLU/SparseLU_column_bmod.h index bf25a33fc..b268c4348 100644 --- a/Eigen/src/SparseLU/SparseLU_column_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_column_bmod.h @@ -46,11 +46,9 @@ * > 0 - number of bytes allocated when run out of space * */ -template -int LU_column_bmod(const int jcol, const int nseg, BlockScalarVector& dense, ScalarVector& tempv, BlockIndexVector& segrep, BlockIndexVector& repfnz, int fpanelc, LU_GlobalLU_t& glu) +template +int SparseLUBase::LU_column_bmod(const int jcol, const int nseg, BlockScalarVector& dense, ScalarVector& tempv, BlockIndexVector& segrep, BlockIndexVector& repfnz, int fpanelc, GlobalLU_t& glu) { - typedef typename IndexVector::Scalar Index; - typedef typename ScalarVector::Scalar Scalar; int jsupno, k, ksub, krep, ksupno; int lptr, nrow, isub, irow, nextlu, new_next, ufirst; int fsupc, nsupc, nsupr, luptr, kfnz, no_zeros; @@ -95,9 +93,6 @@ int LU_column_bmod(const int jcol, const int nseg, BlockScalarVector& dense, Sca nsupr = glu.xlsub(fsupc+1) - glu.xlsub(fsupc); nrow = nsupr - d_fsupc - nsupc; - // NOTE Unlike the original implementation in SuperLU, the only feature - // available here is a sup-col update. - // Perform a triangular solver and block update, // then scatter the result of sup-col update to dense no_zeros = kfnz - fst_col; diff --git a/Eigen/src/SparseLU/SparseLU_column_dfs.h b/Eigen/src/SparseLU/SparseLU_column_dfs.h index 568e0686c..fa8dcf18d 100644 --- a/Eigen/src/SparseLU/SparseLU_column_dfs.h +++ b/Eigen/src/SparseLU/SparseLU_column_dfs.h @@ -42,15 +42,15 @@ * \param m number of rows in the matrix * \param jcol Current column * \param perm_r Row permutation - * \param maxsuper + * \param maxsuper Maximum number of column allowed in a supernode * \param [in,out] nseg Number of segments in current U[*,j] - new segments appended * \param lsub_col defines the rhs vector to start the dfs * \param [in,out] segrep Segment representatives - new segments appended - * \param repfnz + * \param repfnz First nonzero location in each row * \param xprune - * \param marker + * \param marker marker[i] == jj, if i was visited during dfs of current column jj; * \param parent - * \param xplore + * \param xplore working array * \param glu global LU data * \return 0 success * > 0 number of bytes allocated when run out of space @@ -60,6 +60,7 @@ template struct LU_column_dfs_traits { typedef typename IndexVector::Scalar Index; + typedef typename ScalarVector::Scalar Scalar; LU_column_dfs_traits(Index jcol, Index& jsuper, LU_GlobalLU_t& glu) : m_jcol(jcol), m_jsuper_ref(jsuper), m_glu(glu) {} @@ -70,7 +71,7 @@ struct LU_column_dfs_traits void mem_expand(IndexVector& lsub, int& nextl, int chmark) { if (nextl >= m_glu.nzlmax) - LUMemXpand(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions); + SparseLUBase::LUMemXpand(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions); if (chmark != (m_jcol-1)) m_jsuper_ref = IND_EMPTY; } enum { ExpandMem = true }; @@ -80,11 +81,9 @@ struct LU_column_dfs_traits LU_GlobalLU_t& m_glu; }; -template -int LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, int maxsuper, int& nseg, BlockIndexVector& lsub_col, IndexVector& segrep, BlockIndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, LU_GlobalLU_t& glu) +template +int SparseLUBase::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, int maxsuper, int& nseg, BlockIndexVector& lsub_col, IndexVector& segrep, BlockIndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu) { - typedef typename IndexVector::Scalar Index; - typedef typename ScalarVector::Scalar Scalar; int jsuper = glu.supno(jcol); int nextl = glu.xlsub(jcol); diff --git a/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h b/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h index 541785881..d3227469d 100644 --- a/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h +++ b/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h @@ -43,11 +43,9 @@ * > 0 - number of bytes allocated when run out of space * */ -template -int LU_copy_to_ucol(const int jcol, const int nseg, SegRepType& segrep, RepfnzType& repfnz ,IndexVector& perm_r, DenseType& dense, LU_GlobalLU_t& glu) -{ - typedef typename IndexVector::Scalar Index; - typedef typename ScalarVector::Scalar Scalar; +template +int SparseLUBase::LU_copy_to_ucol(const int jcol, const int nseg, IndexVector& segrep, BlockIndexVector& repfnz ,IndexVector& perm_r, BlockScalarVector& dense, GlobalLU_t& glu) +{ Index ksub, krep, ksupno; Index jsupno = glu.supno(jcol); diff --git a/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h b/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h index 1bda70aaf..6d3271aff 100644 --- a/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h +++ b/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h @@ -38,8 +38,8 @@ * \param descendants Number of descendants of each node in the etree * \param relax_end last column in a supernode */ -template -void LU_heap_relax_snode (const int n, IndexVector& et, const int relax_columns, IndexVector& descendants, IndexVector& relax_end) +template +void SparseLUBase::LU_heap_relax_snode (const int n, IndexVector& et, const int relax_columns, IndexVector& descendants, IndexVector& relax_end) { // The etree may not be postordered, but its heap ordered diff --git a/Eigen/src/SparseLU/SparseLU_kernel_bmod.h b/Eigen/src/SparseLU/SparseLU_kernel_bmod.h index d5cad49b1..b15ff9c50 100644 --- a/Eigen/src/SparseLU/SparseLU_kernel_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_kernel_bmod.h @@ -30,7 +30,7 @@ template struct LU_kernel_bmod template EIGEN_DONT_INLINE static void run(const int segsize, BlockScalarVector& dense, ScalarVector& tempv, ScalarVector& lusup, int& luptr, const int nsupr, const int nrow, IndexVector& lsub, const int lptr, const int no_zeros) { - typedef typename ScalarVector::Scalar Scalar; + typedef typename ScalarVector::Scalar Scalar; // First, copy U[*,j] segment from dense(*) to tempv(*) // The result of triangular solve is in tempv[*]; // The result of matric-vector update is in dense[*] diff --git a/Eigen/src/SparseLU/SparseLU_panel_bmod.h b/Eigen/src/SparseLU/SparseLU_panel_bmod.h index 1b31cc31a..ceb6c5938 100644 --- a/Eigen/src/SparseLU/SparseLU_panel_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_panel_bmod.h @@ -34,7 +34,7 @@ * \brief Performs numeric block updates (sup-panel) in topological order. * * Before entering this routine, the original nonzeros in the panel - * were already copied i nto the spa[m,w] ... FIXME to be checked + * were already copied i nto the spa[m,w] * * \param m number of rows in the matrix * \param w Panel size @@ -48,10 +48,9 @@ * * */ -template -void LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, ScalarVector& dense, ScalarVector& tempv, DenseIndexBlock& segrep, DenseIndexBlock& repfnz, LU_perfvalues& perfv, LU_GlobalLU_t& glu) +template +void SparseLUBase::LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, ScalarVector& dense, ScalarVector& tempv, IndexVector& segrep, IndexVector& repfnz, LU_perfvalues& perfv, GlobalLU_t& glu) { - typedef typename ScalarVector::Scalar Scalar; int ksub,jj,nextl_col; int fsupc, nsupc, nsupr, nrow; @@ -190,9 +189,6 @@ void LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, Sca segsize = krep - kfnz + 1; luptr = glu.xlusup(fsupc); - // NOTE : Unlike the original implementation in SuperLU, - // there is no update feature for col-col, 2col-col ... - // Perform a trianglar solve and block update, // then scatter the result of sup-col update to dense[] no_zeros = kfnz - fsupc; @@ -205,4 +201,4 @@ void LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, Sca } // End for each updating supernode } -#endif \ No newline at end of file +#endif diff --git a/Eigen/src/SparseLU/SparseLU_panel_dfs.h b/Eigen/src/SparseLU/SparseLU_panel_dfs.h index 3581f6d9c..164417897 100644 --- a/Eigen/src/SparseLU/SparseLU_panel_dfs.h +++ b/Eigen/src/SparseLU/SparseLU_panel_dfs.h @@ -29,12 +29,12 @@ */ #ifndef SPARSELU_PANEL_DFS_H #define SPARSELU_PANEL_DFS_H - -template -void LU_dfs_kernel(const int jj, IndexVector& perm_r, +template +template +void SparseLUBase::LU_dfs_kernel(const int jj, IndexVector& perm_r, int& nseg, IndexVector& panel_lsub, IndexVector& segrep, - VectorBlock& repfnz_col, IndexVector& xprune, MarkerType& marker, IndexVector& parent, - IndexVector& xplore, LU_GlobalLU_t& glu, + RepfnzType& repfnz_col, IndexVector& xprune, MarkerType& marker, IndexVector& parent, + IndexVector& xplore, GlobalLU_t& glu, int& nextl_col, int krow, Traits& traits ) { @@ -207,8 +207,8 @@ struct LU_panel_dfs_traits Index* m_marker; }; -template -void LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, IndexVector& perm_r, int& nseg, ScalarVector& dense, IndexVector& panel_lsub, IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, LU_GlobalLU_t& glu) +template +void SparseLUBase::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, IndexVector& perm_r, int& nseg, ScalarVector& dense, IndexVector& panel_lsub, IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu) { int nextl_col; // Next available position in panel_lsub[*,jj] diff --git a/Eigen/src/SparseLU/SparseLU_pivotL.h b/Eigen/src/SparseLU/SparseLU_pivotL.h index 4ad49adee..c4a9f1c74 100644 --- a/Eigen/src/SparseLU/SparseLU_pivotL.h +++ b/Eigen/src/SparseLU/SparseLU_pivotL.h @@ -52,12 +52,9 @@ * \return 0 if success, i > 0 if U(i,i) is exactly zero * */ -template -int LU_pivotL(const int jcol, const typename ScalarVector::RealScalar diagpivotthresh, IndexVector& perm_r, IndexVector& iperm_c, int& pivrow, LU_GlobalLU_t& glu) +template +int SparseLUBase::LU_pivotL(const int jcol, const RealScalar diagpivotthresh, IndexVector& perm_r, IndexVector& iperm_c, int& pivrow, GlobalLU_t& glu) { - typedef typename IndexVector::Scalar Index; - typedef typename ScalarVector::Scalar Scalar; - typedef typename ScalarVector::RealScalar RealScalar; Index fsupc = (glu.xsup)((glu.supno)(jcol)); // First column in the supernode containing the column jcol Index nsupc = jcol - fsupc; // Number of columns in the supernode portion, excluding jcol; nsupc >=0 diff --git a/Eigen/src/SparseLU/SparseLU_pruneL.h b/Eigen/src/SparseLU/SparseLU_pruneL.h index f29285bd4..d8c58e039 100644 --- a/Eigen/src/SparseLU/SparseLU_pruneL.h +++ b/Eigen/src/SparseLU/SparseLU_pruneL.h @@ -46,12 +46,9 @@ * \param glu Global LU data * */ -template -void LU_pruneL(const int jcol, const IndexVector& perm_r, const int pivrow, const int nseg, const IndexVector& segrep, BlockIndexVector& repfnz, IndexVector& xprune, LU_GlobalLU_t& glu) +template +void SparseLUBase::LU_pruneL(const int jcol, const IndexVector& perm_r, const int pivrow, const int nseg, const IndexVector& segrep, BlockIndexVector& repfnz, IndexVector& xprune, GlobalLU_t& glu) { - typedef typename IndexVector::Scalar Index; - typedef typename ScalarVector::Scalar Scalar; - // For each supernode-rep irep in U(*,j] int jsupno = glu.supno(jcol); int i,irep,irep1; diff --git a/Eigen/src/SparseLU/SparseLU_relax_snode.h b/Eigen/src/SparseLU/SparseLU_relax_snode.h index a9a0a00c1..8db8619c1 100644 --- a/Eigen/src/SparseLU/SparseLU_relax_snode.h +++ b/Eigen/src/SparseLU/SparseLU_relax_snode.h @@ -37,8 +37,8 @@ * \param descendants Number of descendants of each node in the etree * \param relax_end last column in a supernode */ -template -void LU_relax_snode (const int n, IndexVector& et, const int relax_columns, IndexVector& descendants, IndexVector& relax_end) +template +void SparseLUBase::LU_relax_snode (const int n, IndexVector& et, const int relax_columns, IndexVector& descendants, IndexVector& relax_end) { // compute the number of descendants of each node in the etree diff --git a/Eigen/src/SparseLU/SparseLU_snode_bmod.h b/Eigen/src/SparseLU/SparseLU_snode_bmod.h index 18e6a93d2..beea71e31 100644 --- a/Eigen/src/SparseLU/SparseLU_snode_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_snode_bmod.h @@ -29,11 +29,9 @@ */ #ifndef SPARSELU_SNODE_BMOD_H #define SPARSELU_SNODE_BMOD_H -template -int LU_snode_bmod (const int jcol, const int fsupc, ScalarVector& dense, LU_GlobalLU_t& glu) -{ - typedef typename ScalarVector::Scalar Scalar; - +template +int SparseLUBase::LU_snode_bmod (const int jcol, const int fsupc, ScalarVector& dense, GlobalLU_t& glu) +{ /* lsub : Compressed row subscripts of ( rectangular supernodes ) * xlsub : xlsub[j] is the starting location of the j-th column in lsub(*) * lusup : Numerical values of the rectangular supernodes diff --git a/Eigen/src/SparseLU/SparseLU_snode_dfs.h b/Eigen/src/SparseLU/SparseLU_snode_dfs.h index edb927cdc..199436cd7 100644 --- a/Eigen/src/SparseLU/SparseLU_snode_dfs.h +++ b/Eigen/src/SparseLU/SparseLU_snode_dfs.h @@ -42,55 +42,54 @@ * \param marker (in/out) working vector * \return 0 on success, > 0 size of the memory when memory allocation failed */ - template - int LU_snode_dfs(const int jcol, const int kcol,const MatrixType& mat, IndexVector& xprune, IndexVector& marker, LU_GlobalLU_t& glu) +template +int SparseLUBase::LU_snode_dfs(const int jcol, const int kcol,const MatrixType& mat, IndexVector& xprune, IndexVector& marker, GlobalLU_t& glu) +{ + int mem; + Index nsuper = ++glu.supno(jcol); // Next available supernode number + int nextl = glu.xlsub(jcol); //Index of the starting location of the jcol-th column in lsub + int krow,kmark; + for (int i = jcol; i <=kcol; i++) { - typedef typename IndexVector::Scalar Index; - int mem; - Index nsuper = ++glu.supno(jcol); // Next available supernode number - int nextl = glu.xlsub(jcol); //Index of the starting location of the jcol-th column in lsub - int krow,kmark; - for (int i = jcol; i <=kcol; i++) + // For each nonzero in A(*,i) + for (typename MatrixType::InnerIterator it(mat, i); it; ++it) { - // For each nonzero in A(*,i) - for (typename MatrixType::InnerIterator it(mat, i); it; ++it) + krow = it.row(); + kmark = marker(krow); + if ( kmark != kcol ) { - krow = it.row(); - kmark = marker(krow); - if ( kmark != kcol ) + // First time to visit krow + marker(krow) = kcol; + glu.lsub(nextl++) = krow; + if( nextl >= glu.nzlmax ) { - // First time to visit krow - marker(krow) = kcol; - glu.lsub(nextl++) = krow; - if( nextl >= glu.nzlmax ) - { - mem = LUMemXpand(glu.lsub, glu.nzlmax, nextl, LSUB, glu.num_expansions); - if (mem) return mem; // Memory expansion failed... Return the memory allocated so far - } + mem = LUMemXpand(glu.lsub, glu.nzlmax, nextl, LSUB, glu.num_expansions); + if (mem) return mem; // Memory expansion failed... Return the memory allocated so far } } - glu.supno(i) = nsuper; } - - // If supernode > 1, then make a copy of the subscripts for pruning - if (jcol < kcol) + glu.supno(i) = nsuper; + } + + // If supernode > 1, then make a copy of the subscripts for pruning + if (jcol < kcol) + { + Index new_next = nextl + (nextl - glu.xlsub(jcol)); + while (new_next > glu.nzlmax) { - Index new_next = nextl + (nextl - glu.xlsub(jcol)); - while (new_next > glu.nzlmax) - { - mem = LUMemXpand(glu.lsub, glu.nzlmax, nextl, LSUB, glu.num_expansions); - if (mem) return mem; // Memory expansion failed... Return the memory allocated so far - } - Index ifrom, ito = nextl; - for (ifrom = glu.xlsub(jcol); ifrom < nextl;) - glu.lsub(ito++) = glu.lsub(ifrom++); - for (int i = jcol+1; i <=kcol; i++) glu.xlsub(i) = nextl; - nextl = ito; + mem = LUMemXpand(glu.lsub, glu.nzlmax, nextl, LSUB, glu.num_expansions); + if (mem) return mem; // Memory expansion failed... Return the memory allocated so far } - glu.xsup(nsuper+1) = kcol + 1; // Start of next available supernode - glu.supno(kcol+1) = nsuper; - xprune(kcol) = nextl; - glu.xlsub(kcol+1) = nextl; - return 0; + Index ifrom, ito = nextl; + for (ifrom = glu.xlsub(jcol); ifrom < nextl;) + glu.lsub(ito++) = glu.lsub(ifrom++); + for (int i = jcol+1; i <=kcol; i++) glu.xlsub(i) = nextl; + nextl = ito; } + glu.xsup(nsuper+1) = kcol + 1; // Start of next available supernode + glu.supno(kcol+1) = nsuper; + xprune(kcol) = nextl; + glu.xlsub(kcol+1) = nextl; + return 0; +} #endif \ No newline at end of file -- cgit v1.2.3