SkylineInplaceLU.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Guillaume Saupin <guillaume.saupin@cea.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_SKYLINEINPLACELU_H
00026 #define EIGEN_SKYLINEINPLACELU_H
00027 
00037 template<typename MatrixType>
00038 class SkylineInplaceLU {
00039 protected:
00040     typedef typename MatrixType::Scalar Scalar;
00041     typedef typename MatrixType::Index Index;
00042     
00043     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00044 
00045 public:
00046 
00049     SkylineInplaceLU(MatrixType& matrix, int flags = 0)
00050     : /*m_matrix(matrix.rows(), matrix.cols()),*/ m_flags(flags), m_status(0), m_lu(matrix) {
00051         m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar > ();
00052         m_lu.IsRowMajor ? computeRowMajor() : compute();
00053     }
00054 
00065     void setPrecision(RealScalar v) {
00066         m_precision = v;
00067     }
00068 
00072     RealScalar precision() const {
00073         return m_precision;
00074     }
00075 
00084     void setFlags(int f) {
00085         m_flags = f;
00086     }
00087 
00089     int flags() const {
00090         return m_flags;
00091     }
00092 
00093     void setOrderingMethod(int m) {
00094         m_flags = m;
00095     }
00096 
00097     int orderingMethod() const {
00098         return m_flags;
00099     }
00100 
00102     void compute();
00103     void computeRowMajor();
00104 
00106     //inline const MatrixType& matrixL() const { return m_matrixL; }
00107 
00109     //inline const MatrixType& matrixU() const { return m_matrixU; }
00110 
00111     template<typename BDerived, typename XDerived>
00112     bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x,
00113             const int transposed = 0) const;
00114 
00116     inline bool succeeded(void) const {
00117         return m_succeeded;
00118     }
00119 
00120 protected:
00121     RealScalar m_precision;
00122     int m_flags;
00123     mutable int m_status;
00124     bool m_succeeded;
00125     MatrixType& m_lu;
00126 };
00127 
00131 template<typename MatrixType>
00132 //template<typename _Scalar>
00133 void SkylineInplaceLU<MatrixType>::compute() {
00134     const size_t rows = m_lu.rows();
00135     const size_t cols = m_lu.cols();
00136 
00137     eigen_assert(rows == cols && "We do not (yet) support rectangular LU.");
00138     eigen_assert(!m_lu.IsRowMajor && "LU decomposition does not work with rowMajor Storage");
00139 
00140     for (Index row = 0; row < rows; row++) {
00141         const double pivot = m_lu.coeffDiag(row);
00142 
00143         //Lower matrix Columns update
00144         const Index& col = row;
00145         for (typename MatrixType::InnerLowerIterator lIt(m_lu, col); lIt; ++lIt) {
00146             lIt.valueRef() /= pivot;
00147         }
00148 
00149         //Upper matrix update -> contiguous memory access
00150         typename MatrixType::InnerLowerIterator lIt(m_lu, col);
00151         for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
00152             typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
00153             typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
00154             const double coef = lIt.value();
00155 
00156             uItPivot += (rrow - row - 1);
00157 
00158             //update upper part  -> contiguous memory access
00159             for (++uItPivot; uIt && uItPivot;) {
00160                 uIt.valueRef() -= uItPivot.value() * coef;
00161 
00162                 ++uIt;
00163                 ++uItPivot;
00164             }
00165             ++lIt;
00166         }
00167 
00168         //Upper matrix update -> non contiguous memory access
00169         typename MatrixType::InnerLowerIterator lIt3(m_lu, col);
00170         for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
00171             typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
00172             const double coef = lIt3.value();
00173 
00174             //update lower part ->  non contiguous memory access
00175             for (Index i = 0; i < rrow - row - 1; i++) {
00176                 m_lu.coeffRefLower(rrow, row + i + 1) -= uItPivot.value() * coef;
00177                 ++uItPivot;
00178             }
00179             ++lIt3;
00180         }
00181         //update diag -> contiguous
00182         typename MatrixType::InnerLowerIterator lIt2(m_lu, col);
00183         for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
00184 
00185             typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
00186             typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
00187             const double coef = lIt2.value();
00188 
00189             uItPivot += (rrow - row - 1);
00190             m_lu.coeffRefDiag(rrow) -= uItPivot.value() * coef;
00191             ++lIt2;
00192         }
00193     }
00194 }
00195 
00196 template<typename MatrixType>
00197 void SkylineInplaceLU<MatrixType>::computeRowMajor() {
00198     const size_t rows = m_lu.rows();
00199     const size_t cols = m_lu.cols();
00200 
00201     eigen_assert(rows == cols && "We do not (yet) support rectangular LU.");
00202     eigen_assert(m_lu.IsRowMajor && "You're trying to apply rowMajor decomposition on a ColMajor matrix !");
00203 
00204     for (Index row = 0; row < rows; row++) {
00205         typename MatrixType::InnerLowerIterator llIt(m_lu, row);
00206 
00207 
00208         for (Index col = llIt.col(); col < row; col++) {
00209             if (m_lu.coeffExistLower(row, col)) {
00210                 const double diag = m_lu.coeffDiag(col);
00211 
00212                 typename MatrixType::InnerLowerIterator lIt(m_lu, row);
00213                 typename MatrixType::InnerUpperIterator uIt(m_lu, col);
00214 
00215 
00216                 const Index offset = lIt.col() - uIt.row();
00217 
00218 
00219                 Index stop = offset > 0 ? col - lIt.col() : col - uIt.row();
00220 
00221                 //#define VECTORIZE
00222 #ifdef VECTORIZE
00223                 Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
00224                 Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
00225 
00226 
00227                 Scalar newCoeff = m_lu.coeffLower(row, col) - rowVal.dot(colVal);
00228 #else
00229                 if (offset > 0) //Skip zero value of lIt
00230                     uIt += offset;
00231                 else //Skip zero values of uIt
00232                     lIt += -offset;
00233                 Scalar newCoeff = m_lu.coeffLower(row, col);
00234 
00235                 for (Index k = 0; k < stop; ++k) {
00236                     const Scalar tmp = newCoeff;
00237                     newCoeff = tmp - lIt.value() * uIt.value();
00238                     ++lIt;
00239                     ++uIt;
00240                 }
00241 #endif
00242 
00243                 m_lu.coeffRefLower(row, col) = newCoeff / diag;
00244             }
00245         }
00246 
00247         //Upper matrix update
00248         const Index col = row;
00249         typename MatrixType::InnerUpperIterator uuIt(m_lu, col);
00250         for (Index rrow = uuIt.row(); rrow < col; rrow++) {
00251 
00252             typename MatrixType::InnerLowerIterator lIt(m_lu, rrow);
00253             typename MatrixType::InnerUpperIterator uIt(m_lu, col);
00254             const Index offset = lIt.col() - uIt.row();
00255 
00256             Index stop = offset > 0 ? rrow - lIt.col() : rrow - uIt.row();
00257 
00258 #ifdef VECTORIZE
00259             Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
00260             Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
00261 
00262             Scalar newCoeff = m_lu.coeffUpper(rrow, col) - rowVal.dot(colVal);
00263 #else
00264             if (offset > 0) //Skip zero value of lIt
00265                 uIt += offset;
00266             else //Skip zero values of uIt
00267                 lIt += -offset;
00268             Scalar newCoeff = m_lu.coeffUpper(rrow, col);
00269             for (Index k = 0; k < stop; ++k) {
00270                 const Scalar tmp = newCoeff;
00271                 newCoeff = tmp - lIt.value() * uIt.value();
00272 
00273                 ++lIt;
00274                 ++uIt;
00275             }
00276 #endif
00277             m_lu.coeffRefUpper(rrow, col) = newCoeff;
00278         }
00279 
00280 
00281         //Diag matrix update
00282         typename MatrixType::InnerLowerIterator lIt(m_lu, row);
00283         typename MatrixType::InnerUpperIterator uIt(m_lu, row);
00284 
00285         const Index offset = lIt.col() - uIt.row();
00286 
00287 
00288         Index stop = offset > 0 ? lIt.size() : uIt.size();
00289 #ifdef VECTORIZE
00290         Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
00291         Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
00292         Scalar newCoeff = m_lu.coeffDiag(row) - rowVal.dot(colVal);
00293 #else
00294         if (offset > 0) //Skip zero value of lIt
00295             uIt += offset;
00296         else //Skip zero values of uIt
00297             lIt += -offset;
00298         Scalar newCoeff = m_lu.coeffDiag(row);
00299         for (Index k = 0; k < stop; ++k) {
00300             const Scalar tmp = newCoeff;
00301             newCoeff = tmp - lIt.value() * uIt.value();
00302             ++lIt;
00303             ++uIt;
00304         }
00305 #endif
00306         m_lu.coeffRefDiag(row) = newCoeff;
00307     }
00308 }
00309 
00318 template<typename MatrixType>
00319 template<typename BDerived, typename XDerived>
00320 bool SkylineInplaceLU<MatrixType>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed) const {
00321     const size_t rows = m_lu.rows();
00322     const size_t cols = m_lu.cols();
00323 
00324 
00325     for (Index row = 0; row < rows; row++) {
00326         x->coeffRef(row) = b.coeff(row);
00327         Scalar newVal = x->coeff(row);
00328         typename MatrixType::InnerLowerIterator lIt(m_lu, row);
00329 
00330         Index col = lIt.col();
00331         while (lIt.col() < row) {
00332 
00333             newVal -= x->coeff(col++) * lIt.value();
00334             ++lIt;
00335         }
00336 
00337         x->coeffRef(row) = newVal;
00338     }
00339 
00340 
00341     for (Index col = rows - 1; col > 0; col--) {
00342         x->coeffRef(col) = x->coeff(col) / m_lu.coeffDiag(col);
00343 
00344         const Scalar x_col = x->coeff(col);
00345 
00346         typename MatrixType::InnerUpperIterator uIt(m_lu, col);
00347         uIt += uIt.size()-1;
00348 
00349 
00350         while (uIt) {
00351             x->coeffRef(uIt.row()) -= x_col * uIt.value();
00352             //TODO : introduce --operator
00353             uIt += -1;
00354         }
00355 
00356 
00357     }
00358     x->coeffRef(0) = x->coeff(0) / m_lu.coeffDiag(0);
00359 
00360     return true;
00361 }
00362 
00363 #endif // EIGEN_SKYLINELU_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:23