diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index 0b1347f87..997f4e352 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -153,107 +153,48 @@ class SparseLU for(int j = 0; j < nrhs; ++j) X.col(j) = m_perm_r * B.col(j); - // Forward solve PLy = Pb; - Index fsupc; // First column of the current supernode - Index istart; // Pointer index to the subscript of the current column - Index nsupr; // Number of rows in the current supernode - Index nsupc; // Number of columns in the current supernode - Index nrow; // Number of rows in the non-diagonal part of the supernode - Index luptr; // Pointer index to the current nonzero value - Index iptr; // row index pointer iterator - Index irow; //Current index row - const Scalar * Lval = m_Lstore.valuePtr(); // Nonzero values - Matrix work(n, nrhs); // working vector - work.setZero(); - int j, k, i,jcol; - for (k = 0; k <= m_Lstore.nsuper(); k ++) - { - fsupc = m_Lstore.supToCol()[k]; - istart = m_Lstore.rowIndexPtr()[fsupc]; - nsupr = m_Lstore.rowIndexPtr()[fsupc+1] - istart; - nsupc = m_Lstore.supToCol()[k+1] - fsupc; - nrow = nsupr - nsupc; - - if (nsupc == 1 ) - { - for (j = 0; j < nrhs; j++) - { - luptr = m_Lstore.colIndexPtr()[fsupc]; - for (iptr = istart+1; iptr < m_Lstore.rowIndexPtr()[fsupc+1]; iptr++) - { - irow = m_Lstore.rowIndex()[iptr]; - ++luptr; - X(irow, j) -= X(fsupc, j) * Lval[luptr]; - } - } - } - else - { - // The supernode has more than one column - luptr = m_Lstore.colIndexPtr()[fsupc]; - - // Triangular solve - Map, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); - Map< Matrix, 0, OuterStride<> > U (&(X.data()[fsupc]), nsupc, nrhs, OuterStride<>(n) ); - U = A.template triangularView().solve(U); - - // Matrix-vector product - new (&A) Map, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(nsupr) ); - work.block(0, 0, nrow, nrhs) = A * U; - - //Begin Scatter - for (j = 0; j < nrhs; j++) - { - iptr = istart + nsupc; - for (i = 0; i < nrow; i++) - { - irow = m_Lstore.rowIndex()[iptr]; - X(irow, j) -= work(i, j); // Scatter operation - work(i, j) = Scalar(0); - iptr++; - } - } - } - } // end for all supernodes + //Forward substitution with L + m_Lstore.solveInPlace(X); - // Back solve Ux = y - for (k = m_Lstore.nsuper(); k >= 0; k--) + // Backward solve with U + for (int k = m_Lstore.nsuper(); k >= 0; k--) { - fsupc = m_Lstore.supToCol()[k]; - istart = m_Lstore.rowIndexPtr()[fsupc]; - nsupr = m_Lstore.rowIndexPtr()[fsupc+1] - istart; - nsupc = m_Lstore.supToCol()[k+1] - fsupc; - luptr = m_Lstore.colIndexPtr()[fsupc]; + Index fsupc = m_Lstore.supToCol()[k]; + Index istart = m_Lstore.rowIndexPtr()[fsupc]; + Index nsupr = m_Lstore.rowIndexPtr()[fsupc+1] - istart; + Index nsupc = m_Lstore.supToCol()[k+1] - fsupc; + Index luptr = m_Lstore.colIndexPtr()[fsupc]; if (nsupc == 1) { - for (j = 0; j < nrhs; j++) + for (int j = 0; j < nrhs; j++) { - X(fsupc, j) /= Lval[luptr]; + X(fsupc, j) /= m_Lstore.valuePtr()[luptr]; } } else { - Map, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); - Map< Matrix, 0, OuterStride<> > U (&(X.data()[fsupc]), nsupc, nrhs, OuterStride<>(n) ); + Map, 0, OuterStride<> > A( &(m_Lstore.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); + Map< Matrix, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); U = A.template triangularView().solve(U); } - for (j = 0; j < nrhs; ++j) + for (int j = 0; j < nrhs; ++j) { - for (jcol = fsupc; jcol < fsupc + nsupc; jcol++) + for (int jcol = fsupc; jcol < fsupc + nsupc; jcol++) { - for (i = m_Ustore.outerIndexPtr()[jcol]; i < m_Ustore.outerIndexPtr()[jcol+1]; i++) + typename MappedSparseMatrix::InnerIterator it(m_Ustore, jcol); + for ( ; it; ++it) { - irow = m_Ustore.innerIndexPtr()[i]; - X(irow, j) -= X(jcol, j) * m_Ustore.valuePtr()[i]; + Index irow = it.index(); + X(irow, j) -= X(jcol, j) * it.value(); } } } } // End For U-solve // Permute back the solution - for (j = 0; j < nrhs; ++j) + for (int j = 0; j < nrhs; ++j) X.col(j) = m_perm_c.inverse() * X.col(j); return true; diff --git a/Eigen/src/SparseLU/SparseLU_Matrix.h b/Eigen/src/SparseLU/SparseLU_Matrix.h index 9381189c8..31aeee64d 100644 --- a/Eigen/src/SparseLU/SparseLU_Matrix.h +++ b/Eigen/src/SparseLU/SparseLU_Matrix.h @@ -171,8 +171,13 @@ class SuperNodalMatrix { return m_nsuper; } + class InnerIterator; - class SuperNodeIterator; + template + void solveInPlace( MatrixBase&X) const; + + + protected: Index m_row; // Number of rows @@ -189,7 +194,7 @@ class SuperNodalMatrix }; /** - * \brief InnerIterator class to iterate over nonzero values in the triangular supernodal matrix + * \brief InnerIterator class to iterate over nonzero values of the current column in the supernode * */ template @@ -209,7 +214,7 @@ class SuperNodalMatrix::InnerIterator inline InnerIterator& operator++() { m_idval++; - m_idrow++ ; + m_idrow++; return *this; } inline Scalar value() const { return m_matrix.valuePtr()[m_idval]; } @@ -229,48 +234,80 @@ class SuperNodalMatrix::InnerIterator } 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 + 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 Supernodes in the triangular supernodal matrix + * \brief Solve with the supernode triangular matrix * - * The final goal is to use this class when dealing with supernodes during numerical factorization */ template -class SuperNodalMatrix::SuperNodeIterator +template +void SuperNodalMatrix::solveInPlace( MatrixBase&X) const { - public: - SuperNodeIterator(const SuperNodalMatrix& mat) + Index n = X.rows(); + int nrhs = X.cols(); + const Scalar * Lval = valuePtr(); // Nonzero values + Matrix work(n, nrhs); // working vector + work.setZero(); + for (int k = 0; k <= nsuper(); k ++) { + Index fsupc = supToCol()[k]; // First column of the current supernode + Index istart = rowIndexPtr()[fsupc]; // Pointer index to the subscript of the current column + Index nsupr = rowIndexPtr()[fsupc+1] - istart; // Number of rows in the current supernode + Index nsupc = supToCol()[k+1] - fsupc; // Number of columns in the current supernode + Index nrow = nsupr - nsupc; // Number of rows in the non-diagonal part of the supernode + Index irow; //Current index row - } - 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; - -}; + if (nsupc == 1 ) + { + for (int j = 0; j < nrhs; j++) + { + InnerIterator it(*this, fsupc); + ++it; // Skip the diagonal element + for (; it; ++it) + { + irow = it.row(); + X(irow, j) -= X(fsupc, j) * it.value(); + } + } + } + else + { + // The supernode has more than one column + Index luptr = colIndexPtr()[fsupc]; + + // Triangular solve + Map, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); + Map< Matrix, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); + U = A.template triangularView().solve(U); + + // Matrix-vector product + new (&A) Map, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(nsupr) ); + work.block(0, 0, nrow, nrhs) = A * U; + + //Begin Scatter + for (int j = 0; j < nrhs; j++) + { + Index iptr = istart + nsupc; + for (int i = 0; i < nrow; i++) + { + irow = rowIndex()[iptr]; + X(irow, j) -= work(i, j); // Scatter operation + work(i, j) = Scalar(0); + iptr++; + } + } + } + } +} + + #endif \ No newline at end of file