1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2008 Guillaume Saupin <guillaume.saupin@cea.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_SKYLINEINPLACELU_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_SKYLINEINPLACELU_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup Skyline_Module 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class SkylineInplaceLU 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Inplace LU decomposition of a skyline matrix and associated features 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param MatrixType the type of the matrix of which we are computing the LU factorization 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass SkylineInplaceLU { 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprotected: 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Creates a LU object and compute the respective factorization of \a matrix using 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * flags \a flags. */ 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SkylineInplaceLU(MatrixType& matrix, int flags = 0) 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : /*m_matrix(matrix.rows(), matrix.cols()),*/ m_flags(flags), m_status(0), m_lu(matrix) { 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar > (); 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu.IsRowMajor ? computeRowMajor() : compute(); 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Sets the relative threshold value used to prune zero coefficients during the decomposition. 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Setting a value greater than zero speeds up computation, and yields to an imcomplete 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * factorization with fewer non zero coefficients. Such approximate factors are especially 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * useful to initialize an iterative solver. 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Note that the exact meaning of this parameter might depends on the actual 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * backend. Moreover, not all backends support this feature. 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa precision() */ 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void setPrecision(RealScalar v) { 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_precision = v; 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the current precision. 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa setPrecision() */ 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar precision() const { 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_precision; 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Sets the flags. Possible values are: 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - CompleteFactorization 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - IncompleteFactorization 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - MemoryEfficient 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - one of the ordering methods 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - etc... 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa flags() */ 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void setFlags(int f) { 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_flags = f; 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the current flags */ 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int flags() const { 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_flags; 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void setOrderingMethod(int m) { 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_flags = m; 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int orderingMethod() const { 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_flags; 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Computes/re-computes the LU factorization */ 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute(); 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void computeRowMajor(); 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the lower triangular matrix L */ 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //inline const MatrixType& matrixL() const { return m_matrixL; } 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the upper triangular matrix U */ 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //inline const MatrixType& matrixU() const { return m_matrixU; } 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename BDerived, typename XDerived> 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int transposed = 0) const; 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns true if the factorization succeeded */ 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline bool succeeded(void) const { 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_succeeded; 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprotected: 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar m_precision; 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int m_flags; 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable int m_status; 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_succeeded; 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType& m_lu; 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Computes / recomputes the in place LU decomposition of the SkylineInplaceLU. 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * using the default algorithm. 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath//template<typename _Scalar> 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid SkylineInplaceLU<MatrixType>::compute() { 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const size_t rows = m_lu.rows(); 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const size_t cols = m_lu.cols(); 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rows == cols && "We do not (yet) support rectangular LU."); 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(!m_lu.IsRowMajor && "LU decomposition does not work with rowMajor Storage"); 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index row = 0; row < rows; row++) { 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const double pivot = m_lu.coeffDiag(row); 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //Lower matrix Columns update 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index& col = row; 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (typename MatrixType::InnerLowerIterator lIt(m_lu, col); lIt; ++lIt) { 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lIt.valueRef() /= pivot; 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //Upper matrix update -> contiguous memory access 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerLowerIterator lIt(m_lu, col); 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) { 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerUpperIterator uItPivot(m_lu, row); 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerUpperIterator uIt(m_lu, rrow); 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const double coef = lIt.value(); 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath uItPivot += (rrow - row - 1); 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //update upper part -> contiguous memory access 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (++uItPivot; uIt && uItPivot;) { 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath uIt.valueRef() -= uItPivot.value() * coef; 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++uIt; 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++uItPivot; 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++lIt; 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //Upper matrix update -> non contiguous memory access 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerLowerIterator lIt3(m_lu, col); 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) { 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerUpperIterator uItPivot(m_lu, row); 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const double coef = lIt3.value(); 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //update lower part -> non contiguous memory access 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = 0; i < rrow - row - 1; i++) { 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu.coeffRefLower(rrow, row + i + 1) -= uItPivot.value() * coef; 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++uItPivot; 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++lIt3; 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //update diag -> contiguous 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerLowerIterator lIt2(m_lu, col); 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) { 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerUpperIterator uItPivot(m_lu, row); 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerUpperIterator uIt(m_lu, rrow); 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const double coef = lIt2.value(); 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath uItPivot += (rrow - row - 1); 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu.coeffRefDiag(rrow) -= uItPivot.value() * coef; 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++lIt2; 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid SkylineInplaceLU<MatrixType>::computeRowMajor() { 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const size_t rows = m_lu.rows(); 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const size_t cols = m_lu.cols(); 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rows == cols && "We do not (yet) support rectangular LU."); 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_lu.IsRowMajor && "You're trying to apply rowMajor decomposition on a ColMajor matrix !"); 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index row = 0; row < rows; row++) { 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerLowerIterator llIt(m_lu, row); 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index col = llIt.col(); col < row; col++) { 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_lu.coeffExistLower(row, col)) { 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const double diag = m_lu.coeffDiag(col); 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerLowerIterator lIt(m_lu, row); 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerUpperIterator uIt(m_lu, col); 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index offset = lIt.col() - uIt.row(); 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index stop = offset > 0 ? col - lIt.col() : col - uIt.row(); 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //#define VECTORIZE 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifdef VECTORIZE 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop); 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop); 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar newCoeff = m_lu.coeffLower(row, col) - rowVal.dot(colVal); 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#else 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (offset > 0) //Skip zero value of lIt 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath uIt += offset; 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else //Skip zero values of uIt 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lIt += -offset; 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar newCoeff = m_lu.coeffLower(row, col); 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index k = 0; k < stop; ++k) { 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar tmp = newCoeff; 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath newCoeff = tmp - lIt.value() * uIt.value(); 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++lIt; 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++uIt; 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu.coeffRefLower(row, col) = newCoeff / diag; 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //Upper matrix update 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index col = row; 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerUpperIterator uuIt(m_lu, col); 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index rrow = uuIt.row(); rrow < col; rrow++) { 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerLowerIterator lIt(m_lu, rrow); 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerUpperIterator uIt(m_lu, col); 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index offset = lIt.col() - uIt.row(); 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index stop = offset > 0 ? rrow - lIt.col() : rrow - uIt.row(); 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifdef VECTORIZE 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop); 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop); 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar newCoeff = m_lu.coeffUpper(rrow, col) - rowVal.dot(colVal); 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#else 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (offset > 0) //Skip zero value of lIt 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath uIt += offset; 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else //Skip zero values of uIt 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lIt += -offset; 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar newCoeff = m_lu.coeffUpper(rrow, col); 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index k = 0; k < stop; ++k) { 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar tmp = newCoeff; 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath newCoeff = tmp - lIt.value() * uIt.value(); 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++lIt; 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++uIt; 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu.coeffRefUpper(rrow, col) = newCoeff; 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //Diag matrix update 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerLowerIterator lIt(m_lu, row); 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerUpperIterator uIt(m_lu, row); 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index offset = lIt.col() - uIt.row(); 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index stop = offset > 0 ? lIt.size() : uIt.size(); 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifdef VECTORIZE 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop); 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop); 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar newCoeff = m_lu.coeffDiag(row) - rowVal.dot(colVal); 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#else 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (offset > 0) //Skip zero value of lIt 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath uIt += offset; 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else //Skip zero values of uIt 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lIt += -offset; 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar newCoeff = m_lu.coeffDiag(row); 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index k = 0; k < stop; ++k) { 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar tmp = newCoeff; 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath newCoeff = tmp - lIt.value() * uIt.value(); 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++lIt; 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++uIt; 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu.coeffRefDiag(row) = newCoeff; 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Computes *x = U^-1 L^-1 b 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * If \a transpose is set to SvTranspose or SvAdjoint, the solution 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * of the transposed/adjoint system is computed instead. 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Not all backends implement the solution of the transposed or 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * adjoint system. 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename BDerived, typename XDerived> 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathbool SkylineInplaceLU<MatrixType>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed) const { 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const size_t rows = m_lu.rows(); 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const size_t cols = m_lu.cols(); 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index row = 0; row < rows; row++) { 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x->coeffRef(row) = b.coeff(row); 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar newVal = x->coeff(row); 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerLowerIterator lIt(m_lu, row); 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index col = lIt.col(); 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (lIt.col() < row) { 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath newVal -= x->coeff(col++) * lIt.value(); 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++lIt; 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x->coeffRef(row) = newVal; 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index col = rows - 1; col > 0; col--) { 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x->coeffRef(col) = x->coeff(col) / m_lu.coeffDiag(col); 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar x_col = x->coeff(col); 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::InnerUpperIterator uIt(m_lu, col); 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath uIt += uIt.size()-1; 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (uIt) { 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x->coeffRef(uIt.row()) -= x_col * uIt.value(); 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath //TODO : introduce --operator 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath uIt += -1; 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x->coeffRef(0) = x->coeff(0) / m_lu.coeffDiag(0); 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_SKYLINELU_H 353