1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.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_INCOMPLETE_LUT_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_INCOMPLETE_LUT_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Incomplete LU factorization with dual-threshold strategy 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * During the numerical factorization, two dropping rules are used : 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 1) any element whose magnitude is less than some tolerance is dropped. 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This tolerance is obtained by multiplying the input tolerance @p droptol 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * by the average magnitude of all the original elements in the current row. 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 2) After the elimination of the row, only the @p fill largest elements in 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the L part and the @p fill largest elements in the U part are kept 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * (in addition to the diagonal element ). Note that @p fill is computed from 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the input parameter @p fillfactor which is used the ratio to control the fill_in 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * relatively to the initial number of nonzero elements. 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The two extreme cases are when @p droptol=0 (to keep all the @p fill*2 largest elements) 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and when @p fill=n/2 with @p droptol being different to zero. 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * References : Yousef Saad, ILUT: A dual threshold incomplete LU factorization, 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Numerical Linear Algebra with Applications, 1(4), pp 387-402, 1994. 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * NOTE : The following implementation is derived from the ILUT implementation 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * in the SPARSKIT package, Copyright (C) 2005, the Regents of the University of Minnesota 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * released under the terms of the GNU LGPL: 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * http://www-users.cs.umn.edu/~saad/software/SPARSKIT/README 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * However, Yousef Saad gave us permission to relicense his ILUT code to MPL2. 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * See the Eigen mailing list archive, thread: ILUT, date: July 8, 2012: 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2012/07/msg00064.html 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * alternatively, on GMANE: 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * http://comments.gmane.org/gmane.comp.lib.eigen/3302 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename _Scalar> 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass IncompleteLUT : internal::noncopyable 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _Scalar Scalar; 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,1> Vector; 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef SparseMatrix<Scalar,RowMajor> FactorType; 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef SparseMatrix<Scalar,ColMajor> PermutType; 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename FactorType::Index Index; 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType; 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath IncompleteLUT() 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_droptol(NumTraits<Scalar>::dummy_precision()), m_fillfactor(10), 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false) 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath IncompleteLUT(const MatrixType& mat, RealScalar droptol=NumTraits<Scalar>::dummy_precision(), int fillfactor = 10) 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_droptol(droptol),m_fillfactor(fillfactor), 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_analysisIsOk(false),m_factorizationIsOk(false),m_isInitialized(false) 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(fillfactor != 0); 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(mat); 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows() const { return m_lu.rows(); } 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols() const { return m_lu.cols(); } 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Reports whether previous computation was successful. 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns \c Success if computation was succesful, 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \c NumericalIssue if the matrix.appears to be negative. 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo info() const 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "IncompleteLUT is not initialized."); 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_info; 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void analyzePattern(const MatrixType& amat); 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void factorize(const MatrixType& amat); 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Compute an incomplete LU factorization with dual threshold on the matrix mat 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * No pivoting is done in this version 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath IncompleteLUT<Scalar>& compute(const MatrixType& amat) 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath analyzePattern(amat); 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath factorize(amat); 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_factorizationIsOk == true); 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void setDroptol(RealScalar droptol); 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void setFillfactor(int fillfactor); 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs, typename Dest> 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void _solve(const Rhs& b, Dest& x) const 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x = m_Pinv * b; 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x = m_lu.template triangularView<UnitLower>().solve(x); 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x = m_lu.template triangularView<Upper>().solve(x); 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x = m_P * x; 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs> inline const internal::solve_retval<IncompleteLUT, Rhs> 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solve(const MatrixBase<Rhs>& b) const 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "IncompleteLUT is not initialized."); 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(cols()==b.rows() 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && "IncompleteLUT::solve(): invalid number of rows of the right hand side matrix b"); 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return internal::solve_retval<IncompleteLUT, Rhs>(*this, b.derived()); 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprotected: 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename VectorV, typename VectorI> 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int QuickSplit(VectorV &row, VectorI &ind, int ncut); 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** keeps off-diagonal entries; drops diagonal entries */ 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath struct keep_diag { 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline bool operator() (const Index& row, const Index& col, const Scalar&) const 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return row!=col; 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprotected: 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath FactorType m_lu; 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar m_droptol; 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int m_fillfactor; 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_analysisIsOk; 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_factorizationIsOk; 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo m_info; 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PermutationMatrix<Dynamic,Dynamic,Index> m_P; // Fill-reducing permutation 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // Inverse permutation 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Set control parameter droptol 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param droptol Drop any element whose magnitude is less than this tolerance 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid IncompleteLUT<Scalar>::setDroptol(RealScalar droptol) 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath this->m_droptol = droptol; 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Set control parameter fillfactor 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param fillfactor This is used to compute the number @p fill_in of largest elements to keep on each row. 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid IncompleteLUT<Scalar>::setFillfactor(int fillfactor) 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath this->m_fillfactor = fillfactor; 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Compute a quick-sort split of a vector 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * On output, the vector row is permuted such that its elements satisfy 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * abs(row(i)) >= abs(row(ncut)) if i<ncut 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * abs(row(i)) <= abs(row(ncut)) if i>ncut 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param row The vector of values 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param ind The array of index for the elements in @p row 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param ncut The number of largest elements to keep 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename Scalar> 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename VectorV, typename VectorI> 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathint IncompleteLUT<Scalar>::QuickSplit(VectorV &row, VectorI &ind, int ncut) 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using std::swap; 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int mid; 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int n = row.size(); /* length of the vector */ 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int first, last ; 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ncut--; /* to fit the zero-based indices */ 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath first = 0; 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath last = n-1; 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ncut < first || ncut > last ) return 0; 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath do { 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mid = first; 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar abskey = std::abs(row(mid)); 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int j = first + 1; j <= last; j++) { 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ( std::abs(row(j)) > abskey) { 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++mid; 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath swap(row(mid), row(j)); 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath swap(ind(mid), ind(j)); 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* Interchange for the pivot element */ 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath swap(row(mid), row(first)); 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath swap(ind(mid), ind(first)); 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (mid > ncut) last = mid - 1; 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (mid < ncut ) first = mid + 1; 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } while (mid != ncut ); 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; /* mid is equal to ncut */ 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename Scalar> 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType> 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid IncompleteLUT<Scalar>::analyzePattern(const _MatrixType& amat) 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute the Fill-reducing permutation 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseMatrix<Scalar,ColMajor, Index> mat1 = amat; 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseMatrix<Scalar,ColMajor, Index> mat2 = amat.transpose(); 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Symmetrize the pattern 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // FIXME for a matrix with nearly symmetric pattern, mat2+mat1 is the appropriate choice. 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // on the other hand for a really non-symmetric pattern, mat2*mat1 should be prefered... 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseMatrix<Scalar,ColMajor, Index> AtA = mat2 + mat1; 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath AtA.prune(keep_diag()); 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P); // Then compute the AMD ordering... 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_Pinv = m_P.inverse(); // ... and the inverse permutation 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_analysisIsOk = true; 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename Scalar> 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType> 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using std::sqrt; 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using std::swap; 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using std::abs; 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert((amat.rows() == amat.cols()) && "The factorization should be done on a square matrix"); 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int n = amat.cols(); // Size of the matrix 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu.resize(n,n); 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Declare Working vectors and variables 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Vector u(n) ; // real values of the row -- maximum size is n -- 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXi ju(n); // column position of the values in u -- maximum size is n 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXi jr(n); // Indicate the position of the nonzero elements in the vector u -- A zero location is indicated by -1 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Apply the fill-reducing permutation 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SparseMatrix<Scalar,RowMajor, Index> mat; 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat = amat.twistedBy(m_Pinv); 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Initialization 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jr.fill(-1); 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ju.fill(0); 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u.fill(0); 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // number of largest elements to keep in each row: 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int fill_in = static_cast<int> (amat.nonZeros()*m_fillfactor)/n+1; 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (fill_in > n) fill_in = n; 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // number of largest nonzero elements to keep in the L and the U part of the current row: 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int nnzL = fill_in/2; 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int nnzU = nnzL; 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu.reserve(n * (nnzL + nnzU + 1)); 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // global loop over the rows of the sparse matrix 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int ii = 0; ii < n; ii++) 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // 1 - copy the lower and the upper part of the row i of mat in the working vector u 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int sizeu = 1; // number of nonzero elements in the upper part of the current row 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int sizel = 0; // number of nonzero elements in the lower part of the current row 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ju(ii) = ii; 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u(ii) = 0; 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jr(ii) = ii; 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar rownorm = 0; 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename FactorType::InnerIterator j_it(mat, ii); // Iterate through the current row ii 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (; j_it; ++j_it) 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int k = j_it.index(); 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k < ii) 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // copy the lower part 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ju(sizel) = k; 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u(sizel) = j_it.value(); 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jr(k) = sizel; 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++sizel; 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (k == ii) 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u(ii) = j_it.value(); 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // copy the upper part 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int jpos = ii + sizeu; 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ju(jpos) = k; 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u(jpos) = j_it.value(); 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jr(k) = jpos; 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++sizeu; 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rownorm += internal::abs2(j_it.value()); 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // 2 - detect possible zero row 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(rownorm==0) 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = NumericalIssue; 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return; 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Take the 2-norm of the current row as a relative tolerance 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rownorm = sqrt(rownorm); 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // 3 - eliminate the previous nonzero rows 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int jj = 0; 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int len = 0; 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (jj < sizel) 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // In order to eliminate in the correct order, 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // we must select first the smallest column index among ju(jj:sizel) 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int k; 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int minrow = ju.segment(jj,sizel-jj).minCoeff(&k); // k is relative to the segment 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath k += jj; 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (minrow != ju(jj)) 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // swap the two locations 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int j = ju(jj); 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath swap(ju(jj), ju(k)); 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jr(minrow) = jj; jr(j) = k; 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath swap(u(jj), u(k)); 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Reset this location 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jr(minrow) = -1; 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Start elimination 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename FactorType::InnerIterator ki_it(m_lu, minrow); 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (ki_it && ki_it.index() < minrow) ++ki_it; 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_internal_assert(ki_it && ki_it.col()==minrow); 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar fact = u(jj) / ki_it.value(); 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // drop too small elements 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(abs(fact) <= m_droptol) 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jj++; 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath continue; 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // linear combination of the current row ii and the row minrow 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++ki_it; 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (; ki_it; ++ki_it) 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar prod = fact * ki_it.value(); 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int j = ki_it.index(); 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int jpos = jr(j); 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (jpos == -1) // fill-in element 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int newpos; 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (j >= ii) // dealing with the upper part 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath newpos = ii + sizeu; 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sizeu++; 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_internal_assert(sizeu<=n); 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else // dealing with the lower part 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath newpos = sizel; 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sizel++; 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_internal_assert(sizel<=ii); 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ju(newpos) = j; 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u(newpos) = -prod; 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jr(j) = newpos; 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u(jpos) -= prod; 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // store the pivot element 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u(len) = fact; 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ju(len) = minrow; 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++len; 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jj++; 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } // end of the elimination on the row ii 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // reset the upper part of the pointer jr to zero 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1; 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // 4 - partially sort and insert the elements in the m_lu matrix 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // sort the L-part of the row 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sizel = len; 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath len = (std::min)(sizel, nnzL); 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Vector::SegmentReturnType ul(u.segment(0, sizel)); 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename VectorXi::SegmentReturnType jul(ju.segment(0, sizel)); 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath QuickSplit(ul, jul, len); 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // store the largest m_fill elements of the L part 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu.startVec(ii); 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int k = 0; k < len; k++) 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k); 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // store the diagonal element 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // apply a shifting rule to avoid zero pivots (we are doing an incomplete factorization) 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (u(ii) == Scalar(0)) 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u(ii) = sqrt(m_droptol) * rownorm; 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu.insertBackByOuterInnerUnordered(ii, ii) = u(ii); 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // sort the U-part of the row 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // apply the dropping rule first 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath len = 0; 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int k = 1; k < sizeu; k++) 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(abs(u(ii+k)) > m_droptol * rownorm ) 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++len; 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u(ii + len) = u(ii + k); 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ju(ii + len) = ju(ii + k); 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sizeu = len + 1; // +1 to take into account the diagonal element 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath len = (std::min)(sizeu, nnzU); 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1)); 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename VectorXi::SegmentReturnType juu(ju.segment(ii+1, sizeu-1)); 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath QuickSplit(uu, juu, len); 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // store the largest elements of the U part 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int k = ii + 1; k < ii + len; k++) 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k); 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu.finalize(); 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu.makeCompressed(); 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_factorizationIsOk = true; 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = Success; 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, typename Rhs> 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct solve_retval<IncompleteLUT<_MatrixType>, Rhs> 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : solve_retval_base<IncompleteLUT<_MatrixType>, Rhs> 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef IncompleteLUT<_MatrixType> Dec; 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Dest> void evalTo(Dest& dst) const 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dec()._solve(rhs(),dst); 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_INCOMPLETE_LUT_H 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 467