diff options
Diffstat (limited to 'Eigen/src')
-rw-r--r-- | Eigen/src/SparseLU/SparseLU.h | 126 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_Structs.h | 10 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_column_bmod.h | 58 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_column_dfs.h | 30 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_copy_to_ucol.h | 25 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_panel_bmod.h | 46 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_panel_dfs.h | 43 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_pivotL.h | 22 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_pruneL.h | 28 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_snode_bmod.h | 20 |
10 files changed, 219 insertions, 189 deletions
diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index a4b4fa98b..36b1ce570 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -27,19 +27,12 @@ #define EIGEN_SPARSE_LU namespace Eigen { - -template <typename _MatrixType> -class SparseLU; -#include <Ordering.h> + +// Data structure needed by all routines #include <SparseLU_Structs.h> -#include <SparseLU_Memory.h> -#include <SparseLU_Utils.h> -#include <SuperNodalMatrix.h> +#include <SparseLU_Matrix.h> -#include <SparseLU_Coletree.h> -#include <SparseLU_heap_relax_snode.h> -#include <SparseLU_relax_snode.h> /** * \ingroup SparseLU_Module * \brief Sparse supernodal LU factorization for general matrices @@ -62,7 +55,7 @@ class SparseLU typedef Matrix<Index,Dynamic,1> IndexVector; typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; public: - SparseLU():m_isInitialized(true),m_symmetricmode(false),m_fact(DOFACT),m_diagpivotthresh(1.0) + SparseLU():m_isInitialized(true),m_symmetricmode(false),m_diagpivotthresh(1.0) { initperfvalues(); } @@ -106,7 +99,7 @@ class SparseLU } - /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. + /** \returns the solution X of \f$ A X = b \f$ using the current decomposition of A. * * \sa compute() */ @@ -122,20 +115,34 @@ class SparseLU protected: // Functions void initperfvalues(); - template <typename IndexVector, typename ScalarVector> - int LU_snode_dfs(const int jcol, const int kcol, const IndexVector* asub, - const IndexVector* colptr, IndexVector& xprune, IndexVector& marker, LU_GlobalLU_t& glu); - - template <typename Index, typename ScalarVector> + int LU_snode_dfs(const int jcol, const int kcol, const IndexVector* asub, + const IndexVector* colptr, IndexVector& xprune, IndexVector& marker, LU_GlobalLU_t& glu); int LU_dsnode_bmod (const Index jcol, const Index jsupno, const Index fsupc, - ScalarVector& dense, ScalarVector& tempv, LU_GlobalLu_t& Glu); + ScalarVector& dense, LU_GlobalLU_t& Glu); + int LU_pivotL(const int jcol, const RealScalar diagpivotthresh, IndexVector& perm_r, + IndexVector& iperm_c, int& pivrow, GlobalLU_t& Glu); + 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); + void LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, + ScalarVector& dense, ScalarVector& tempv, IndexVector& segrep, + IndexVector& repfnz, LU_GlobalLU_t& glu); + int LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, IndexVector& nseg, + IndexVector& lsub_col, IndexVector& segrep, IndexVector& repfnz, + IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, LU_GlobalLU_t& glu); + int LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense, ScalarVector& tempv, + IndexVector& segrep, IndexVector& repfnz, int fpanelc, LU_GlobalLU_t& Glu); + int LU_copy_to_ucol(const int jcol, const int nseg, IndexVector& segrep, IndexVector& repfnz, + IndexVector& perm_r, ScalarVector& dense, LU_GlobalLU_t& glu); + void LU_pruneL(const int jcol, const IndexVector& perm_r, const int pivrow, const int nseg, + const IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, GlobalLU_t& Glu) // Variables mutable ComputationInfo m_info; bool m_isInitialized; bool m_factorizationIsOk; bool m_analysisIsOk; - fact_t m_fact; NCMatrix m_mat; // The input (permuted ) matrix SCMatrix m_Lstore; // The lower triangular matrix (supernodal) NCMatrix m_Ustore; // The upper triangular matrix @@ -146,7 +153,8 @@ class SparseLU ScalarVector m_work; // Scalar work vector IndexVector m_iwork; //Index work vector static LU_GlobalLU_t m_glu; // persistent data to facilitate multiple factors - // should be defined as a class member + // FIXME All fields of this struct can be defined separately as class members + // SuperLU/SparseLU options bool m_symmetricmode; @@ -179,7 +187,10 @@ void SparseLU::initperfvalues() m_fillfactor = 20; } - +// Functions needed by the anaysis phase +#include <SparseLU_Coletree.h> +// Ordering interface +#include <Ordering.h> /** * Compute the column permutation to minimize the fill-in (file amd.c ) * @@ -206,7 +217,7 @@ void SparseLU::analyzePattern(const MatrixType& mat) // Apply the permutation to the column of the input matrix - m_mat = mat * m_perm_c; //FIXME Check if this is valid, check as well how to permute only the index + m_mat = mat * m_perm_c; // Compute the column elimination tree of the permuted matrix if (m_etree.size() == 0) m_etree.resize(m_mat.cols()); @@ -234,6 +245,21 @@ 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_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> + /** * - Numerical factorization * - Interleaved with the symbolic factorization @@ -284,7 +310,7 @@ void SparseLU::factorize(const MatrixType& matrix) idx += m; VectorBlock<IndexVector> xplore(m_iwork, idx, m); idx += m; - VectorBlock<IndexVector> repnfnz(m_iwork, idx, maxpanel); + VectorBlock<IndexVector> repfnz(m_iwork, idx, maxpanel); idx += maxpanel; VectorBlock<IndexVector> panel_lsub(m_iwork, idx, maxpanel) idx += maxpanel; @@ -324,8 +350,6 @@ void SparseLU::factorize(const MatrixType& matrix) supno(0) = IND_EMPTY; xsup(0) = xlsub(0) = xusub(0) = xlusup(0) = 0; - int panel_size = m_panel_size; - int wdef = m_panel_size; // upper bound on panel width // 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 @@ -348,9 +372,9 @@ void SparseLU::factorize(const MatrixType& matrix) info = LU_snode_dfs(jcol, kcol, m_mat.innerIndexPtr(), m_mat.outerIndexPtr(), xprune, marker); if ( info ) { + std::cerr << "MEMORY ALLOCATION FAILED IN SNODE_DFS() \n"; m_info = NumericalIssue; m_factorizationIsOk = false; - std::cerr << "MEMORY ALLOCATION FAILED IN SNODE_DFS() \n"; return; } nextu = xusub(jcol); //starting location of column jcol in ucol @@ -377,13 +401,14 @@ void SparseLU::factorize(const MatrixType& matrix) dense(it.row()) = it.val(); // Numeric update within the snode - LU_snode_bmod(icol, jsupno, fsupc, dense, tempv); + LU_snode_bmod(icol, jsupno, fsupc, dense, glu); // Eliminate the current column info = LU_pivotL(icol, m_diagpivotthresh, m_perm_r, m_iperm_c, pivrow, m_glu); if ( info ) { m_info = NumericalIssue; + std::cerr<< "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT " << info <<std::endl; m_factorizationIsOk = false; return; } @@ -394,7 +419,7 @@ void SparseLU::factorize(const MatrixType& matrix) { // Work on one panel of panel_size columns // Adjust panel size so that a panel won't overlap with the next relaxed snode. - panel_size = w_def; + int panel_size = wdef; // upper bound on panel width for (k = jcol + 1; k < std::min(jcol+panel_size, n); k++) { if (relax_end(k) != IND_EMPTY) @@ -419,31 +444,33 @@ void SparseLU::factorize(const MatrixType& matrix) nseg = nseg1; // begin after all the panel segments //Depth-first-search for the current column - VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m); //FIXME - VectorBlock<IndexVector> repfnz_k(repfnz, k, m); //FIXME + VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m); + VectorBlock<IndexVector> repfnz_k(repfnz, k, m); info = LU_column_dfs(m, jj, perm_r, nseg, panel_lsub(k), segrep, repfnz_k, xprune, marker, parent, xplore, m_glu); if ( !info ) { + std::cerr << "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() \n"; m_info = NumericalIssue; m_factorizationIsOk = false; return; } // Numeric updates to this column - VectorBlock<IndexVector> dense_k(dense, k, m); //FIXME - VectorBlock<IndexVector> segrep_k(segrep, nseg1, m) // FIXME Check the length + VectorBlock<IndexVector> dense_k(dense, k, m); + VectorBlock<IndexVector> segrep_k(segrep, nseg1, m); info = 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"; m_info = NumericalIssue; m_factorizationIsOk = false; return; } // Copy the U-segments to ucol(*) - //FIXME Check that repfnz_k, dense_k... have stored references to modified columns info = LU_copy_to_col(jj, nseg, segrep, repfnz_k, perm_r, dense_k, m_glu); if ( info ) { + std::cerr << "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() \n"; m_info = NumericalIssue; m_factorizationIsOk = false; return; @@ -453,6 +480,7 @@ void SparseLU::factorize(const MatrixType& matrix) info = LU_pivotL(jj, m_diagpivotthresh, m_perm_r, iperm_c, pivrow, m_glu); if ( info ) { + std::cerr<< "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT " << info <<std::endl; m_info = NumericalIssue; m_factorizationIsOk = false; return; @@ -472,7 +500,7 @@ void SparseLU::factorize(const MatrixType& matrix) } // end else } // end for -- end elimination - // Adjust row permutation in the case of rectangular matrices + // Adjust row permutation in the case of rectangular matrices... Deprecated if (m > n ) { k = 0; @@ -504,18 +532,18 @@ void SparseLU::factorize(const MatrixType& matrix) } template<typename Rhs, typename Dest> -bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const +bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &X) const { eigen_assert(m_isInitialized && "The matrix should be factorized first"); EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0, THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); - x = b; /* on return, x is overwritten by the computed solution */ + X = b; /* on return, X is overwritten by the computed solution */ int nrhs = b.cols(); // Permute the right hand side to form Pr*B - x = m_perm_r * x; + X = m_perm_r * X; // Forward solve PLy = Pb; Index fsupc; // First column of the current supernode @@ -547,7 +575,7 @@ bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const { irow = m_Lstore.rowIndex()[iptr]; ++luptr; - x(irow, j) -= x(fsupc, j) * Lval[luptr]; + X(irow, j) -= X(fsupc, j) * Lval[luptr]; } } } @@ -558,8 +586,8 @@ bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const // Triangular solve luptr = m_Lstore.colIndexPtr()[fsupc]; //FIXME Should be outside the loop Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); -// Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride > u( &(x(fsupc,0)), nsupc, nrhs, OuterStride<>(x.rows()) ); - Matrix<Scalar,Dynamic,Dynamic>& u = x.block(fsupc, 0, nsupc, nrhs); //FIXME Check this +// Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride > u( &(X(fsupc,0)), nsupc, nrhs, OuterStride<>(X.rows()) ); + Matrix<Scalar,Dynamic,Dynamic>& u = X.block(fsupc, 0, nsupc, nrhs); //FIXME Check this u = A.triangularView<Lower>().solve(u); // Matrix-vector product @@ -573,7 +601,7 @@ bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const for (i = 0; i < nrow; i++) { irow = m_Lstore.rowIndex()[iptr]; - x(irow, j) -= work(i, j); // Scatter operation + X(irow, j) -= work(i, j); // Scatter operation work(i, j) = Scalar(0); iptr++; } @@ -594,13 +622,13 @@ bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const { for (j = 0; j < nrhs; j++) { - x(fsupc, j) /= Lval[luptr]; + X(fsupc, j) /= Lval[luptr]; } } else { Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); - Matrix<Scalar,Dynamic,Dynamic>& u = x.block(fsupc, 0, nsupc, nrhs); + Matrix<Scalar,Dynamic,Dynamic>& u = X.block(fsupc, 0, nsupc, nrhs); u = A.triangularView<Upper>().solve(u); } @@ -608,17 +636,17 @@ bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const { for (jcol = fsupc; jcol < fsupc + nsupc; jcol++) { - for (i = m_Ustore.outerIndexPtr()[jcol]; i < m_Ustore.outerIndexPtr()[jcol]; i++) - { - irow = m_Ustore.InnerIndices()[i]; - x(irow, j) -= x(irow, jcol) * m_Ustore.Values()[i]; - } + for (i = m_Ustore.outerIndexPtr()[jcol]; i < m_Ustore.outerIndexPtr()[jcol]; i++) + { + irow = m_Ustore.InnerIndices()[i]; + X(irow, j) -= X(irow, jcol) * m_Ustore.Values()[i]; + } } } } // End For U-solve // Permute back the solution - x = x * m_perm_c; + X = m_perm_c * X; return true; } diff --git a/Eigen/src/SparseLU/SparseLU_Structs.h b/Eigen/src/SparseLU/SparseLU_Structs.h index 1394eccdf..618d05eac 100644 --- a/Eigen/src/SparseLU/SparseLU_Structs.h +++ b/Eigen/src/SparseLU/SparseLU_Structs.h @@ -82,19 +82,12 @@ */ #ifndef EIGEN_LU_STRUCTS #define EIGEN_LU_STRUCTS -namespace Eigen { - -#define LU_NBR_MEMTYPE 4 /* 0: lusup - 1: ucol - 2: lsub - 3: usub */ -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; template <typename ScalarVector, typename IndexVector> struct { + typedef typename IndexVector::Index Index; 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 @@ -113,5 +106,4 @@ struct { int num_expansions; } GlobalLU_t; -}// End namespace Eigen #endif
\ No newline at end of file diff --git a/Eigen/src/SparseLU/SparseLU_column_bmod.h b/Eigen/src/SparseLU/SparseLU_column_bmod.h index bed4f9519..965a0c0ad 100644 --- a/Eigen/src/SparseLU/SparseLU_column_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_column_bmod.h @@ -54,35 +54,40 @@ * \param segrep segment representative ... * \param repfnz ??? First nonzero column in each row ??? ... * \param fpanelc First column in the current panel - * \param Glu Global LU data. + * \param glu Global LU data. * \return 0 - successful return * > 0 - number of bytes allocated when run out of space * */ template <typename ScalarVector, typename IndexVector> -int SparseLU::LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense, ScalarVector& tempv, IndexVector& segrep, IndexVector& repfnz, int fpanelc, LU_GlobalLu_t& Glu) +int SparseLU::LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense, ScalarVector& tempv, IndexVector& segrep, IndexVector& repfnz, int fpanelc, LU_GlobalLU_t& glu) { - int jsupno, k, ksub, krep, krep_ind, ksupno; + int jsupno, k, ksub, krep, krep_ind, ksupno; + int fsupc, nsupc, nsupr, luptr, kfnz, no_zeros; /* krep = representative of current k-th supernode * fsupc = first supernodal column * nsupc = number of columns in a supernode * nsupr = number of rows in a supernode * luptr = location of supernodal LU-block in storage * kfnz = first nonz in the k-th supernodal segment - * no-zeros = no lf leading zeros in a supernodal U-segment + * no_zeros = no lf leading zeros in a supernodal U-segment */ - IndexVector& xsup = Glu.xsup; - IndexVector& supno = Glu.supno; - IndexVector& lsub = Glu.lsub; - IndexVector& xlsub = Glu.xlsub; - IndexVector& xlusup = Glu.xlusup; - ScalarVector& lusup = Glu.lusup; - Index& nzlumax = Glu.nzlumax; + IndexVector& xsup = glu.xsup; + IndexVector& supno = glu.supno; + IndexVector& lsub = glu.lsub; + IndexVector& xlsub = glu.xlsub; + IndexVector& xlusup = glu.xlusup; + ScalarVector& lusup = glu.lusup; + Index& nzlumax = glu.nzlumax; int jsupno = supno(jcol); // For each nonzero supernode segment of U[*,j] in topological order k = nseg - 1; + int d_fsupc; // distance between the first column of the current panel and the + // first column of the current snode + int fst_col; // First column within small LU update + int segsize; for (ksub = 0; ksub < nseg; ksub++) { krep = segrep(k); k--; @@ -110,35 +115,36 @@ int SparseLU::LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense krep_ind = lptr + nsupc - 1; // NOTE Unlike the original implementation in SuperLU, the only feature - // here is a sup-col update. + // 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; // First, copy U[*,j] segment from dense(*) to tempv(*) isub = lptr + no_zeros; - for (i = 0; i ww segsize; i++) + for (i = 0; i < segsize; i++) { irow = lsub(isub); - tempv(i) = densee(irow); + tempv(i) = dense(irow); ++isub; } // Dense triangular solve -- start effective triangle luptr += nsupr * no_zeros + no_zeros; // Form Eigen matrix and vector Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) ); - Map<ScalarVector> u(tempv.data(), segsize); + VectorBlock<ScalarVector> u(tempv, 0, segsize); + u = A.triangularView<Lower>().solve(u); // Dense matrix-vector product y <-- A*x luptr += segsize; - new (&A) (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(nsupr) ); - Map<ScalarVector> l( &(tempv.data()[segsize]), segsize); + new (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(nsupr) ); + VectorBlock<ScalarVector> l(tempv, segsize, nrow); l= A * u; // Scatter tempv[] into SPA dense[] as a temporary storage isub = lptr + no_zeros; - for (i = 0; i w segsize; i++) + for (i = 0; i < segsize; i++) { irow = lsub(isub); dense(irow) = tempv(i); @@ -150,8 +156,8 @@ int SparseLU::LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense for (i = 0; i < nrow; i++) { irow = lsub(isub); - dense(irow) -= tempv(segsize + i); - tempv(segsize + i) = Scalar(0.0); + dense(irow) -= l(i); + l(i) = Scalar(0.0); ++isub; } } // end if jsupno @@ -165,9 +171,9 @@ int SparseLU::LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense new_next = nextlu + xlsub(fsupc + 1) - xlsub(fsupc); while (new_next > nzlumax ) { - mem = LUmemXpand<Scalar>(Glu.lusup, nzlumax, nextlu, LUSUP, Glu); + mem = LUmemXpand<Scalar>(glu.lusup, nzlumax, nextlu, LUSUP, glu); if (mem) return mem; - lsub = Glu.lsub; //FIXME Why is it updated here. + //lsub = glu.lsub; // Should not be updated here } for (isub = xlsub(fsupc); isub < xlsub(fsupc+1); isub++) @@ -183,8 +189,8 @@ int SparseLU::LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense /* For more updates within the panel (also within the current supernode), * should start from the first column of the panel, or the first column * of the supernode, whichever is bigger. There are two cases: - * 1) fsupc < fpanelc, then fst_col <- fpanelc - * 2) fsupc >= fpanelc, then fst_col <-fsupc + * 1) fsupc < fpanelc, then fst_col <-- fpanelc + * 2) fsupc >= fpanelc, then fst_col <-- fsupc */ fst_col = std::max(fsupc, fpanelc); @@ -203,11 +209,11 @@ int SparseLU::LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense // points to the beginning of jcol in snode L\U(jsupno) ufirst = xlusup(jcol) + d_fsupc; Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); - Map<ScalarVector> l( &(lusup.data()[ufirst]), nsupc ); + VectorBlock<ScalarVector> u(lusup, ufirst, nsupc); u = A.triangularView().solve(u); new (&A) Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr+nsupc]), nrow, nsupc, OuterStride<>(nsupr) ); - Map<ScalarVector> l( &(lusup.data()[ufirst+nsupc]), nsupr ); + VectorBlock<ScalarVector> l(lusup, ufirst+nsupc, nrow); l = l - A * u; } // End if fst_col diff --git a/Eigen/src/SparseLU/SparseLU_column_dfs.h b/Eigen/src/SparseLU/SparseLU_column_dfs.h index 1c832d60e..7fda536a9 100644 --- a/Eigen/src/SparseLU/SparseLU_column_dfs.h +++ b/Eigen/src/SparseLU/SparseLU_column_dfs.h @@ -65,13 +65,13 @@ * \param marker * \param parent * \param xplore - * \param Glu global LU data + * \param glu global LU data * \return 0 success * > 0 number of bytes allocated when run out of space * */ -template <typename IndexVector> -int SparseLU::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, IndexVector& nseg IndexVector& lsub_col, IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, LU_GlobalLu_t& Glu) +template <typename IndexVector, typename ScalarVector> +int SparseLU::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, IndexVector& nseg IndexVector& lsub_col, IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, LU_GlobalLU_t& glu) { typedef typename IndexVector::IndexVector; @@ -85,15 +85,15 @@ int SparseLU::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, In int xdfs, maxdfs, kpar; int mem; // Initialize pointers - IndexVector& xsup = Glu.xsup; - IndexVector& supno = Glu.supno; - IndexVector& lsub = Glu.lsub; - IndexVector& xlsub = Glu.xlsub; - IndexVector& nzlmax = Glu.nzlmax; + IndexVector& xsup = glu.xsup; + IndexVector& supno = glu.supno; + IndexVector& lsub = glu.lsub; + IndexVector& xlsub = glu.xlsub; + IndexVector& nzlmax = glu.nzlmax; nsuper = supno(jcol); jsuper = nsuper; - nextl = xlsup(jcol); + nextl = xlsub(jcol); VectorBlock<IndexVector> marker2(marker, 2*m, m); // For each nonzero in A(*,jcol) do dfs for (k = 0; lsub_col[k] != IND_EMPTY; k++) @@ -106,16 +106,16 @@ int SparseLU::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, In if (kmark == jcol) continue; // For each unmarker nbr krow of jcol - // krow is in L: place it in structure of L(*,jcol) marker2(krow) = jcol; kperm = perm_r(krow); if (kperm == IND_EMPTY ) { + // krow is in L: place it in structure of L(*,jcol) lsub(nextl++) = krow; // krow is indexed into A if ( nextl >= nzlmax ) { - mem = LUMemXpand<IndexVector>(lsub, nzlmax, nextl, LSUB, Glu); + mem = LUMemXpand<IndexVector>(lsub, nzlmax, nextl, LSUB, glu); if ( mem ) return mem; } if (kmark != jcolm1) jsuper = IND_EMPTY; // Row index subset testing @@ -157,13 +157,13 @@ int SparseLU::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, In marker2(kchild) = jcol; chperm = perm_r(kchild); - // if kchild is in L: place it in L(*,k) if (chperm == IND_EMPTY) { + // if kchild is in L: place it in L(*,k) lsub(nextl++) = kchild; if (nextl >= nzlmax) { - mem = LUMemXpand<IndexVector>(lsub, nzlmax, nextl, LSUB); + mem = LUMemXpand<IndexVector>(lsub, nzlmax, nextl, LSUB, glu); if (mem) return mem; } if (chmark != jcolm1) jsuper = IND_EMPTY; @@ -201,7 +201,7 @@ int SparseLU::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, In // place supernode-rep krep in postorder DFS. // backtrack dfs to its parent - segrep(nseg) = ;krep; + segrep(nseg) = krep; ++nseg; kpar = parent(krep); // Pop from stack, mimic recursion if (kpar == IND_EMPTY) break; // dfs done @@ -217,7 +217,7 @@ int SparseLU::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, In } // for each nonzero ... - // check to see if j belongs in the same supeprnode as j-1 + // check to see if j belongs in the same supernode as j-1 if ( jcol == 0 ) { // Do nothing for column 0 nsuper = supno(0) = 0 ; diff --git a/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h b/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h index dc53edcfb..c97bc6aa4 100644 --- a/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h +++ b/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h @@ -50,28 +50,28 @@ * \param jcol current column to update * \param nseg Number of segments in the U part * \param segrep segment representative ... - * \param repfnz ??? First nonzero column in each row ??? ... + * \param repfnz First nonzero column in each row ... * \param perm_r Row permutation * \param dense Store the full representation of the column - * \param Glu Global LU data. + * \param glu Global LU data. * \return 0 - successful return * > 0 - number of bytes allocated when run out of space * */ template <typename ScalarVector, typename IndexVector> -int SparseLU::LU_copy_to_ucol(const int jcol, const int nseg, IndexVector& segrep, IndexVector& repfnz, IndexVector& perm_r, ScalarVector& dense, LU_GlobalLu_t& Glu) +int SparseLU::LU_copy_to_ucol(const int jcol, const int nseg, IndexVector& segrep, IndexVector& repfnz, IndexVector& perm_r, ScalarVector& dense, LU_GlobalLU_t& glu) { Index ksupno, k, ksub, krep, ksupno; typedef typename IndexVector::Index; - IndexVector& xsup = Glu.xsup; - IndexVector& supno = Glu.supno; - IndexVector& lsub = Glu.lsub; - IndexVector& xlsub = Glu.xlsub; + IndexVector& xsup = glu.xsup; + IndexVector& supno = glu.supno; + IndexVector& lsub = glu.lsub; + IndexVector& xlsub = glu.xlsub; ScalarVector& ucol = GLu.ucol; - IndexVector& usub = Glu.usub; - IndexVector& xusub = Glu.xusub; - Index& nzumax = Glu.nzumax; + IndexVector& usub = glu.usub; + IndexVector& xusub = glu.xusub; + Index& nzumax = glu.nzumax; Index jsupno = supno(jcol); @@ -95,12 +95,11 @@ int SparseLU::LU_copy_to_ucol(const int jcol, const int nseg, IndexVector& segre new_next = nextu + segsize; while (new_next > nzumax) { - mem = LU_MemXpand<ScalarVector>(ucol, nzumax, nextu, UCOL, Glu); + mem = LU_MemXpand<ScalarVector>(ucol, nzumax, nextu, UCOL, glu); if (mem) return mem; - mem = LU_MemXpand<Index>(usub, nzumax, nextu, USUB, Glu); + mem = LU_MemXpand<Index>(usub, nzumax, nextu, USUB, glu); if (mem) return mem; - lsub = Glu.lsub; //FIXME Why setting this as well ?? } for (i = 0; i < segsize; i++) diff --git a/Eigen/src/SparseLU/SparseLU_panel_bmod.h b/Eigen/src/SparseLU/SparseLU_panel_bmod.h index 93daa938c..212ecfa6a 100644 --- a/Eigen/src/SparseLU/SparseLU_panel_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_panel_bmod.h @@ -56,21 +56,21 @@ * \param nseg Number of segments in the U part * \param dense Store the full representation of the panel * \param tempv working array - * \param segrep in ... - * \param repfnz in ... - * \param Glu Global LU data. + * \param segrep segment representative... first row in the segment + * \param repfnz First nonzero rows + * \param glu Global LU data. * * */ -template <typename VectorType> -void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, VectorType& dense, VectorType& tempv, VectorXi& segrep, VectorXi& repfnz, LU_GlobalLu_t& Glu) +template <typename IndexVector, typename ScalarVector> +void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, ScalarVector& dense, ScalarVector& tempv, IndexVector& segrep, IndexVector& repfnz, LU_GlobalLU_t& glu) { - VectorXi& xsup = Glu.xsup; - VectorXi& supno = Glu.supno; - VectorXi& lsub = Glu.lsub; - VectorXi& xlsub = Glu.xlsub; - VectorXi& xlusup = Glu.xlusup; - VectorType& lusup = Glu.lusup; + IndexVector& xsup = glu.xsup; + IndexVector& supno = glu.supno; + IndexVector& lsub = glu.lsub; + IndexVector& xlsub = glu.xlsub; + IndexVector& xlusup = glu.xlusup; + ScalarVector& lusup = glu.lusup; int i,ksub,jj,nextl_col,irow; int fsupc, nsupc, nsupr, nrow; @@ -96,10 +96,7 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int nrow = nsupr - nsupc; lptr = xlsub(fsupc); krep_ind = lptr + nsupc - 1; - - repfnz_col = repfnz; - dense_col = dense; - + // NOTE : Unlike the original implementation in SuperLU, the present implementation // does not include a 2-D block update. @@ -107,8 +104,8 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int for (jj = jcol; jj < jcol + w; jj++) { nextl_col = (jj-jcol) * m; - VectorBlock<VectorXi> repfnz_col(repfnz.segment(nextl_col, m)); // First nonzero column index for each row - VectorBLock<VectorXi> dense_col(dense.segment(nextl_col, m)); // Scatter/gather entire matrix column from/to here + VectorBlock<IndexVector> repfnz_col(repfnz.segment(nextl_col, m)); // First nonzero column index for each row + VectorBLock<IndexVector> dense_col(dense.segment(nextl_col, m)); // Scatter/gather entire matrix column from/to here kfnz = repfnz_col(krep); if ( kfnz == IND_EMPTY ) @@ -123,8 +120,7 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int // Perform a trianglar solve and block update, // then scatter the result of sup-col update to dense[] no_zeros = kfnz - fsupc; - - // Copy U[*,j] segment from dense[*] to tempv[*] : + // 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_col[*] isub = lptr + no_zeros; @@ -138,19 +134,21 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int luptr += nsupr * no_zeros + no_zeros; // triangular solve with Eigen Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) ); - Map<Matrix<Scalar,Dynamic,1> > u( tempv.data(), segsize); +// Map<Matrix<Scalar,Dynamic,1> > u( tempv.data(), segsize); + VectorBlock<ScalarVector> u(tempv, 0, segsize); u = A.triangularView<Lower>().solve(u); luptr += segsize; // Dense Matrix vector product y <-- A*x; new (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(nsupr) ); - Map<VectorType> l( &(tempv.data()[segsize]), segsize); +// Map<ScalarVector> l( &(tempv.data()[segsize]), nrow); + VectorBlock<ScalarVector> l(tempv, segsize, nrow); l= A * u; // Scatter tempv(*) into SPA dense(*) such that tempv(*) // can be used for the triangular solve of the next // column of the panel. The y will be copied into ucol(*) - // after the whole panel has been finished. + // after the whole panel has been finished... after column_dfs() and column_bmod() isub = lptr + no_zeros; for (i = 0; i < segsize; i++) @@ -166,8 +164,8 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int for (i = 0; i < nrow; i++) { irow = lsub(isub); - dense_col(irow) -= tempv(segsize + i); - tempv(segsize + i) = 0; + dense_col(irow) -= l(i); + l(i) = Scalar(0); ++isub; } diff --git a/Eigen/src/SparseLU/SparseLU_panel_dfs.h b/Eigen/src/SparseLU/SparseLU_panel_dfs.h index 97e5121db..d3c2906b2 100644 --- a/Eigen/src/SparseLU/SparseLU_panel_dfs.h +++ b/Eigen/src/SparseLU/SparseLU_panel_dfs.h @@ -62,45 +62,50 @@ * marker[i] == jj, if i was visited during dfs of current column jj; * marker1[i] >= jcol, if i was visited by earlier columns in this panel; * - * \param m number of rows in the matrix - * \param w Panel size - * \param jcol Starting column of the panel - * \param A Input matrix in column-major storage - * \param perm_r Row permutation - * \param nseg Number of U segments - * ... + * \param [in]m number of rows in the matrix + * \param [in]w Panel size + * \param [in]jcol Starting column of the panel + * \param [in]A Input matrix in column-major storage + * \param [in]perm_r Row permutation + * \param [out]nseg Number of U segments + * \param [out]dense Accumulate the column vectors of the panel + * \param [out]panel_lsub Subscripts of the row in the panel + * \param [out]segrep Segment representative i.e first nonzero row of each segment + * \param [out]repfnz First nonzero location in each row + * \param [out]xprune + * \param [out]marker + * * */ -template <typename MatrixType, typename VectorType> -void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, VectorXi& perm_r, int& nseg, VectorType& dense, VectorXi& panel_lsub, VectorXi& segrep, VectorXi& repfnz, VectorXi& xprune, VectorXi& marker, VectorXi& parent, VectorXi& xplore, LU_GlobalLu_t& Glu) +template <typename MatrixType, typename IndexVector, typename ScalarVector> +void SparseLU::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) { int jj; // Index through each column in the panel int nextl_col; // Next available position in panel_lsub[*,jj] int krow; // Row index of the current element int kperm; // permuted row index - int krep; // Supernode reprentative of the current row + int krep; // Supernode representative of the current row int kmark; int chperm, chmark, chrep, oldrep, kchild; int myfnz; // First nonzero element in the current column int xdfs, maxdfs, kpar; // Initialize pointers -// VectorXi& marker1 = marker.block(m, m); - VectorBlock<VectorXi> marker1(marker, m, m); +// IndexVector& marker1 = marker.block(m, m); + VectorBlock<IndexVector> marker1(marker, m, m); nseg = 0; - VectorXi& xsup = Glu.xsup; - VectorXi& supno = Glu.supno; - VectorXi& lsub = Glu.lsub; - VectorXi& xlsub = Glu.xlsub; + IndexVector& xsup = Glu.xsup; + IndexVector& supno = Glu.supno; + IndexVector& lsub = Glu.lsub; + IndexVector& xlsub = Glu.xlsub; // For each column in the panel for (jj = jcol; jj < jcol + w; jj++) { nextl_col = (jj - jcol) * m; - //FIXME - VectorBlock<VectorXi> repfnz_col(repfnz.segment(nextl_col, m)); // First nonzero location in each row - VectorBlock<VectorXi> dense_col(dense.segment(nextl_col, m)); // Accumulate a column vector here + VectorBlock<IndexVector> repfnz_col(repfnz, nextl_col, m); // First nonzero location in each row + VectorBlock<IndexVector> dense_col(dense,nextl_col, m); // Accumulate a column vector here // For each nnz in A[*, jj] do depth first search diff --git a/Eigen/src/SparseLU/SparseLU_pivotL.h b/Eigen/src/SparseLU/SparseLU_pivotL.h index 3bfe14e7e..32da92481 100644 --- a/Eigen/src/SparseLU/SparseLU_pivotL.h +++ b/Eigen/src/SparseLU/SparseLU_pivotL.h @@ -67,28 +67,30 @@ * \return 0 if success, i > 0 if U(i,i) is exactly zero * */ -template <typename Scalar> -int SparseLU::LU_pivotL(const int jcol, const Scalar u, VectorXi& perm_r, VectorXi& iperm_c, int& pivrow, GlobalLU_t& Glu) +template <typename IndexVector, typename ScalarVector> +int SparseLU::LU_pivotL(const int jcol, const RealScalar diagpivotthresh, IndexVector& perm_r, IndexVector& iperm_c, int& pivrow, GlobalLU_t& Glu) { + typedef typename IndexVector::Index Index; + typedef typename ScalarVector::Scalar Scalar; // Initialize pointers - VectorXi& lsub = Glu.lsub; // Compressed row subscripts of ( rectangular supernodes ??) - VectorXi& xlsub = Glu.xlsub; // xlsub[j] is the starting location of the j-th column in lsub(*) - Scalar* lusup = Glu.lusup.data(); // Numerical values of the rectangular supernodes - VectorXi& xlusup = Glu.xlusup; // xlusup[j] is the starting location of the j-th column in lusup(*) + IndexVector& lsub = Glu.lsub; // Compressed row subscripts of L rectangular supernodes. + IndexVector& xlsub = Glu.xlsub; // pointers to the beginning of each column subscript in lsub + ScalarVector& lusup = Glu.lusup; // Numerical values of L ordered by columns + IndexVector& xlusup = Glu.xlusup; // pointers to the beginning of each colum in lusup 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 Index lptr = xlsub(fsupc); // pointer to the starting location of the row subscripts for this supernode portion Index nsupr = xlsub(fsupc+1) - lptr; // Number of rows in the supernode - Scalar* lu_sup_ptr = &(lusup[xlusup(fsupc)]); // Start of the current supernode - Scalar* lu_col_ptr = &(lusup[xlusup(jcol)]); // Start of jcol in the supernode + Scalar* lu_sup_ptr = &(lusup.data()[xlusup(fsupc)]); // Start of the current supernode + Scalar* lu_col_ptr = &(lusup.data()[xlusup(jcol)]); // Start of jcol in the supernode Index* lsub_ptr = &(lsub.data()[lptr]); // Start of row indices of the supernode // Determine the largest abs numerical value for partial pivoting Index diagind = iperm_c(jcol); // diagonal index Scalar pivmax = 0.0; Index pivptr = nsupc; - Index diag = -1; + Index diag = IND_EMPTY; Index old_pivptr = nsupc; Scalar rtemp; for (isub = nsupc; isub < nsupr; ++isub) { @@ -127,7 +129,7 @@ int SparseLU::LU_pivotL(const int jcol, const Scalar u, VectorXi& perm_r, Vector // Interchange row subscripts if (pivptr != nsupc ) { - std::swap( lsub_ptr(pivptr), lsub_ptr(nsupc) ); + std::swap( lsub_ptr[pivptr], lsub_ptr[nsupc] ); // Interchange numerical values as well, for the two rows in the whole snode // such that L is indexed the same way as A for (icol = 0; icol <= nsupc; icol++) diff --git a/Eigen/src/SparseLU/SparseLU_pruneL.h b/Eigen/src/SparseLU/SparseLU_pruneL.h index 687717d52..dd092b778 100644 --- a/Eigen/src/SparseLU/SparseLU_pruneL.h +++ b/Eigen/src/SparseLU/SparseLU_pruneL.h @@ -47,35 +47,35 @@ /** * \brief Prunes the L-structure. * - * It prunes the L-structure of supernodes whose L-structure constains the current pivot row "pivrow" + * It prunes the L-structure of supernodes whose L-structure contains the current pivot row "pivrow" * * * \param jcol The current column of L * \param [in]perm_r Row permutation * \param [out]pivrow The pivot row - * \param nseg Number of segments ??? + * \param nseg Number of segments * \param segrep * \param repfnz * \param [out]xprune * \param Glu Global LU data * */ -template <typename VectorType> -void SparseLU::LU_pruneL(const int jcol, const VectorXi& perm_r, const int pivrow, const int nseg, const VectorXi& segrep, VectorXi& repfnz, VectorXi& xprune, GlobalLU_t& Glu) +template <typename IndexVector, typename ScalarVector> +void SparseLU::LU_pruneL(const int jcol, const IndexVector& perm_r, const int pivrow, const int nseg, const IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, GlobalLU_t& Glu) { // Initialize pointers - VectorXi& xsup = Glu.xsup; - VectorXi& supno = Glu.supno; - VectorXi& lsub = Glu.lsub; - VectorXi& xlsub = Glu.xlsub; - VectorType& lusup = Glu.lusup; - VectorXi& xlusup = Glu.xlusup; + IndexVector& xsup = Glu.xsup; + IndexVector& supno = Glu.supno; + IndexVector& lsub = Glu.lsub; + IndexVector& xlsub = Glu.xlsub; + ScalarVector& lusup = Glu.lusup; + IndexVector& xlusup = Glu.xlusup; // For each supernode-rep irep in U(*,j] int jsupno = supno(jcol); int i,irep,irep1; bool movnum, do_prune = false; - int kmin, kmax, ktemp, minloc, maxloc; + Index kmin, kmax, ktemp, minloc, maxloc; for (i = 0; i < nseg; i++) { irep = segrep(i); @@ -125,9 +125,7 @@ void SparseLU::LU_pruneL(const int jcol, const VectorXi& perm_r, const int pivro { // kmin below pivrow (not yet pivoted), and kmax // above pivrow: interchange the two suscripts - ktemp = lsub(kmin); - lsub(kmin) = lsub(kmax); - lsub(kmax) = ktemp; + std::swap(lsub(kmin), lsub(kmax)); // If the supernode has only one column, then we // only keep one set of subscripts. For any subscript @@ -144,7 +142,7 @@ void SparseLU::LU_pruneL(const int jcol, const VectorXi& perm_r, const int pivro } } // end while - xprune(irep) = kmin; + xprune(irep) = kmin; //Pruning } // end if do_prune } // end pruning } // End for each U-segment diff --git a/Eigen/src/SparseLU/SparseLU_snode_bmod.h b/Eigen/src/SparseLU/SparseLU_snode_bmod.h index 6130a5622..1d6bed8bb 100644 --- a/Eigen/src/SparseLU/SparseLU_snode_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_snode_bmod.h @@ -47,13 +47,13 @@ namespace internal { #define SPARSELU_SNODE_BMOD_H template <typename Index, typename ScalarVector> int SparseLU::LU_dsnode_bmod (const Index jcol, const Index jsupno, const Index fsupc, - ScalarVector& dense, ScalarVector& tempv, LU_GlobalLu_t& Glu) + ScalarVector& dense, LU_GlobalLU_t& glu) { typedef typename Matrix<Index, Dynamic, Dynamic> IndexVector; - IndexVector& lsub = Glu.lsub; // Compressed row subscripts of ( rectangular supernodes ??) - IndexVector& xlsub = Glu.xlsub; // xlsub[j] is the starting location of the j-th column in lsub(*) - ScalarVector& lusup = Glu.lusup; // Numerical values of the rectangular supernodes - IndexVector& xlusup = Glu.xlusup; // xlusup[j] is the starting location of the j-th column in lusup(*) + IndexVector& lsub = glu.lsub; // Compressed row subscripts of ( rectangular supernodes ??) + IndexVector& xlsub = glu.xlsub; // xlsub[j] is the starting location of the j-th column in lsub(*) + ScalarVector& lusup = glu.lusup; // Numerical values of the rectangular supernodes + IndexVector& xlusup = glu.xlusup; // xlusup[j] is the starting location of the j-th column in lusup(*) int nextlu = xlusup(jcol); // Starting location of the next column to add int irow, isub; @@ -75,14 +75,16 @@ int SparseLU::LU_dsnode_bmod (const Index jcol, const Index jsupno, const Index int nrow = nsupr - nsupc; // Number of rows in the off-diagonal blocks - // Solve the triangular system for U(fsupc:jcol, jcol) with L(fspuc..., fsupc:jcol) + // Solve the triangular system for U(fsupc:jcol, jcol) with L(fspuc:jcol, fsupc:jcol) Map<Matrix<Scalar,Dynamic,Dynamic>,0,OuterStride<> > A( &(lusup.data()[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); - Map<Matrix<Scalar,Dynamic,1> > u(&(lusup.data()[ufirst]), nsupc); - u = A.triangularView<Lower>().solve(u); +// Map<Matrix<Scalar,Dynamic,1> > u(&(lusup.data()[ufirst]), nsupc); + VectorBlock<ScalarVector> u(lusup, ufirst, nsupc); + u = A.triangularView<Lower>().solve(u); // Call the Eigen dense triangular solve interface // 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], nsupc); +// Map<Matrix<Scalar,Dynamic,1> > l(&(lusup.data()[ufirst+nsupc], nrow); + VectorBlock<ScalarVector> l(lusup, ufirst+nsupc, nrow); l = l - A * u; return 0; |