src/Skyline/SkylineInplaceLU.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008 Guillaume Saupin <guillaume.saupin@cea.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_SKYLINEINPLACELU_H
0011 #define EIGEN_SKYLINEINPLACELU_H
0012
0013 namespace Eigen {
0014
0015 /** \ingroup Skyline_Module
0016  *
0017  * \class SkylineInplaceLU
0018  *
0019  * \brief Inplace LU decomposition of a skyline matrix and associated features
0020  *
0021  * \param MatrixType the type of the matrix of which we are computing the LU factorization
0022  *
0023  */
0024 template<typename MatrixType>
0025 class SkylineInplaceLU {
0026 protected:
0027     typedef typename MatrixType::Scalar Scalar;
0028     typedef typename MatrixType::Index Index;
0029
0030     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
0031
0032 public:
0033
0034     /** Creates a LU object and compute the respective factorization of \a matrix using
0035      * flags \a flags. */
0036     SkylineInplaceLU(MatrixType& matrix, int flags = 0)
0037     : /*m_matrix(matrix.rows(), matrix.cols()),*/ m_flags(flags), m_status(0), m_lu(matrix) {
0038         m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar > ();
0039         m_lu.IsRowMajor ? computeRowMajor() : compute();
0040     }
0041
0042     /** Sets the relative threshold value used to prune zero coefficients during the decomposition.
0043      *
0044      * Setting a value greater than zero speeds up computation, and yields to an incomplete
0045      * factorization with fewer non zero coefficients. Such approximate factors are especially
0046      * useful to initialize an iterative solver.
0047      *
0048      * Note that the exact meaning of this parameter might depends on the actual
0049      * backend. Moreover, not all backends support this feature.
0050      *
0051      * \sa precision() */
0052     void setPrecision(RealScalar v) {
0053         m_precision = v;
0054     }
0055
0056     /** \returns the current precision.
0057      *
0058      * \sa setPrecision() */
0059     RealScalar precision() const {
0060         return m_precision;
0061     }
0062
0063     /** Sets the flags. Possible values are:
0064      *  - CompleteFactorization
0065      *  - IncompleteFactorization
0066      *  - MemoryEfficient
0067      *  - one of the ordering methods
0068      *  - etc...
0069      *
0070      * \sa flags() */
0071     void setFlags(int f) {
0072         m_flags = f;
0073     }
0074
0075     /** \returns the current flags */
0076     int flags() const {
0077         return m_flags;
0078     }
0079
0080     void setOrderingMethod(int m) {
0081         m_flags = m;
0082     }
0083
0084     int orderingMethod() const {
0085         return m_flags;
0086     }
0087
0088     /** Computes/re-computes the LU factorization */
0089     void compute();
0090     void computeRowMajor();
0091
0092     /** \returns the lower triangular matrix L */
0093     //inline const MatrixType& matrixL() const { return m_matrixL; }
0094
0095     /** \returns the upper triangular matrix U */
0096     //inline const MatrixType& matrixU() const { return m_matrixU; }
0097
0098     template<typename BDerived, typename XDerived>
0099     bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x,
0100             const int transposed = 0) const;
0101
0102     /** \returns true if the factorization succeeded */
0103     inline bool succeeded(void) const {
0104         return m_succeeded;
0105     }
0106
0107 protected:
0108     RealScalar m_precision;
0109     int m_flags;
0110     mutable int m_status;
0111     bool m_succeeded;
0112     MatrixType& m_lu;
0113 };
0114
0115 /** Computes / recomputes the in place LU decomposition of the SkylineInplaceLU.
0116  * using the default algorithm.
0117  */
0118 template<typename MatrixType>
0119 //template<typename _Scalar>
0120 void SkylineInplaceLU<MatrixType>::compute() {
0121     const size_t rows = m_lu.rows();
0122     const size_t cols = m_lu.cols();
0123
0124     eigen_assert(rows == cols && "We do not (yet) support rectangular LU.");
0125     eigen_assert(!m_lu.IsRowMajor && "LU decomposition does not work with rowMajor Storage");
0126
0127     for (Index row = 0; row < rows; row++) {
0128         const double pivot = m_lu.coeffDiag(row);
0129
0130         //Lower matrix Columns update
0131         const Index& col = row;
0132         for (typename MatrixType::InnerLowerIterator lIt(m_lu, col); lIt; ++lIt) {
0133             lIt.valueRef() /= pivot;
0134         }
0135
0136         //Upper matrix update -> contiguous memory access
0137         typename MatrixType::InnerLowerIterator lIt(m_lu, col);
0138         for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
0139             typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
0140             typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
0141             const double coef = lIt.value();
0142
0143             uItPivot += (rrow - row - 1);
0144
0145             //update upper part  -> contiguous memory access
0146             for (++uItPivot; uIt && uItPivot;) {
0147                 uIt.valueRef() -= uItPivot.value() * coef;
0148
0149                 ++uIt;
0150                 ++uItPivot;
0151             }
0152             ++lIt;
0153         }
0154
0155         //Upper matrix update -> non contiguous memory access
0156         typename MatrixType::InnerLowerIterator lIt3(m_lu, col);
0157         for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
0158             typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
0159             const double coef = lIt3.value();
0160
0161             //update lower part ->  non contiguous memory access
0162             for (Index i = 0; i < rrow - row - 1; i++) {
0163                 m_lu.coeffRefLower(rrow, row + i + 1) -= uItPivot.value() * coef;
0164                 ++uItPivot;
0165             }
0166             ++lIt3;
0167         }
0168         //update diag -> contiguous
0169         typename MatrixType::InnerLowerIterator lIt2(m_lu, col);
0170         for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
0171
0172             typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
0173             typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
0174             const double coef = lIt2.value();
0175
0176             uItPivot += (rrow - row - 1);
0177             m_lu.coeffRefDiag(rrow) -= uItPivot.value() * coef;
0178             ++lIt2;
0179         }
0180     }
0181 }
0182
0183 template<typename MatrixType>
0184 void SkylineInplaceLU<MatrixType>::computeRowMajor() {
0185     const size_t rows = m_lu.rows();
0186     const size_t cols = m_lu.cols();
0187
0188     eigen_assert(rows == cols && "We do not (yet) support rectangular LU.");
0189     eigen_assert(m_lu.IsRowMajor && "You're trying to apply rowMajor decomposition on a ColMajor matrix !");
0190
0191     for (Index row = 0; row < rows; row++) {
0192         typename MatrixType::InnerLowerIterator llIt(m_lu, row);
0193
0194
0195         for (Index col = llIt.col(); col < row; col++) {
0196             if (m_lu.coeffExistLower(row, col)) {
0197                 const double diag = m_lu.coeffDiag(col);
0198
0199                 typename MatrixType::InnerLowerIterator lIt(m_lu, row);
0200                 typename MatrixType::InnerUpperIterator uIt(m_lu, col);
0201
0202
0203                 const Index offset = lIt.col() - uIt.row();
0204
0205
0206                 Index stop = offset > 0 ? col - lIt.col() : col - uIt.row();
0207
0208                 //#define VECTORIZE
0209 #ifdef VECTORIZE
0210                 Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
0211                 Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
0212
0213
0214                 Scalar newCoeff = m_lu.coeffLower(row, col) - rowVal.dot(colVal);
0215 #else
0216                 if (offset > 0) //Skip zero value of lIt
0217                     uIt += offset;
0218                 else //Skip zero values of uIt
0219                     lIt += -offset;
0220                 Scalar newCoeff = m_lu.coeffLower(row, col);
0221
0222                 for (Index k = 0; k < stop; ++k) {
0223                     const Scalar tmp = newCoeff;
0224                     newCoeff = tmp - lIt.value() * uIt.value();
0225                     ++lIt;
0226                     ++uIt;
0227                 }
0228 #endif
0229
0230                 m_lu.coeffRefLower(row, col) = newCoeff / diag;
0231             }
0232         }
0233
0234         //Upper matrix update
0235         const Index col = row;
0236         typename MatrixType::InnerUpperIterator uuIt(m_lu, col);
0237         for (Index rrow = uuIt.row(); rrow < col; rrow++) {
0238
0239             typename MatrixType::InnerLowerIterator lIt(m_lu, rrow);
0240             typename MatrixType::InnerUpperIterator uIt(m_lu, col);
0241             const Index offset = lIt.col() - uIt.row();
0242
0243             Index stop = offset > 0 ? rrow - lIt.col() : rrow - uIt.row();
0244
0245 #ifdef VECTORIZE
0246             Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
0247             Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
0248
0249             Scalar newCoeff = m_lu.coeffUpper(rrow, col) - rowVal.dot(colVal);
0250 #else
0251             if (offset > 0) //Skip zero value of lIt
0252                 uIt += offset;
0253             else //Skip zero values of uIt
0254                 lIt += -offset;
0255             Scalar newCoeff = m_lu.coeffUpper(rrow, col);
0256             for (Index k = 0; k < stop; ++k) {
0257                 const Scalar tmp = newCoeff;
0258                 newCoeff = tmp - lIt.value() * uIt.value();
0259
0260                 ++lIt;
0261                 ++uIt;
0262             }
0263 #endif
0264             m_lu.coeffRefUpper(rrow, col) = newCoeff;
0265         }
0266
0267
0268         //Diag matrix update
0269         typename MatrixType::InnerLowerIterator lIt(m_lu, row);
0270         typename MatrixType::InnerUpperIterator uIt(m_lu, row);
0271
0272         const Index offset = lIt.col() - uIt.row();
0273
0274
0275         Index stop = offset > 0 ? lIt.size() : uIt.size();
0276 #ifdef VECTORIZE
0277         Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
0278         Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
0279         Scalar newCoeff = m_lu.coeffDiag(row) - rowVal.dot(colVal);
0280 #else
0281         if (offset > 0) //Skip zero value of lIt
0282             uIt += offset;
0283         else //Skip zero values of uIt
0284             lIt += -offset;
0285         Scalar newCoeff = m_lu.coeffDiag(row);
0286         for (Index k = 0; k < stop; ++k) {
0287             const Scalar tmp = newCoeff;
0288             newCoeff = tmp - lIt.value() * uIt.value();
0289             ++lIt;
0290             ++uIt;
0291         }
0292 #endif
0293         m_lu.coeffRefDiag(row) = newCoeff;
0294     }
0295 }
0296
0297 /** Computes *x = U^-1 L^-1 b
0298  *
0299  * If \a transpose is set to SvTranspose or SvAdjoint, the solution
0300  * of the transposed/adjoint system is computed instead.
0301  *
0302  * Not all backends implement the solution of the transposed or
0303  * adjoint system.
0304  */
0305 template<typename MatrixType>
0306 template<typename BDerived, typename XDerived>
0307 bool SkylineInplaceLU<MatrixType>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed) const {
0308     const size_t rows = m_lu.rows();
0309     const size_t cols = m_lu.cols();
0310
0311
0312     for (Index row = 0; row < rows; row++) {
0313         x->coeffRef(row) = b.coeff(row);
0314         Scalar newVal = x->coeff(row);
0315         typename MatrixType::InnerLowerIterator lIt(m_lu, row);
0316
0317         Index col = lIt.col();
0318         while (lIt.col() < row) {
0319
0320             newVal -= x->coeff(col++) * lIt.value();
0321             ++lIt;
0322         }
0323
0324         x->coeffRef(row) = newVal;
0325     }
0326
0327
0328     for (Index col = rows - 1; col > 0; col--) {
0329         x->coeffRef(col) = x->coeff(col) / m_lu.coeffDiag(col);
0330
0331         const Scalar x_col = x->coeff(col);
0332
0333         typename MatrixType::InnerUpperIterator uIt(m_lu, col);
0334         uIt += uIt.size()-1;
0335
0336
0337         while (uIt) {
0338             x->coeffRef(uIt.row()) -= x_col * uIt.value();
0339             //TODO : introduce --operator
0340             uIt += -1;
0341         }
0342
0343
0344     }
0345     x->coeffRef(0) = x->coeff(0) / m_lu.coeffDiag(0);
0346
0347     return true;
0348 }
0349
0350 } // end namespace Eigen
0351
0352 #endif // EIGEN_SKYLINEINPLACELU_H