src/SparseCore/TriangularSolver.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.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_SPARSETRIANGULARSOLVER_H
0011 #define EIGEN_SPARSETRIANGULARSOLVER_H
0012
0013 namespace Eigen {
0014
0015 namespace internal {
0016
0017 template<typename Lhs, typename Rhs, int Mode,
0018   int UpLo = (Mode & Lower)
0019            ? Lower
0020            : (Mode & Upper)
0021            ? Upper
0022            : -1,
0023   int StorageOrder = int(traits<Lhs>::Flags) & RowMajorBit>
0024 struct sparse_solve_triangular_selector;
0025
0026 // forward substitution, row-major
0027 template<typename Lhs, typename Rhs, int Mode>
0028 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,RowMajor>
0029 {
0030   typedef typename Rhs::Scalar Scalar;
0031   typedef evaluator<Lhs> LhsEval;
0032   typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
0033   static void run(const Lhs& lhs, Rhs& other)
0034   {
0035     LhsEval lhsEval(lhs);
0036     for(Index col=0 ; col<other.cols() ; ++col)
0037     {
0038       for(Index i=0; i<lhs.rows(); ++i)
0039       {
0040         Scalar tmp = other.coeff(i,col);
0041         Scalar lastVal(0);
0042         Index lastIndex = 0;
0043         for(LhsIterator it(lhsEval, i); it; ++it)
0044         {
0045           lastVal = it.value();
0046           lastIndex = it.index();
0047           if(lastIndex==i)
0048             break;
0049           tmp -= lastVal * other.coeff(lastIndex,col);
0050         }
0051         if (Mode & UnitDiag)
0052           other.coeffRef(i,col) = tmp;
0053         else
0054         {
0055           eigen_assert(lastIndex==i);
0056           other.coeffRef(i,col) = tmp/lastVal;
0057         }
0058       }
0059     }
0060   }
0061 };
0062
0063 // backward substitution, row-major
0064 template<typename Lhs, typename Rhs, int Mode>
0065 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,RowMajor>
0066 {
0067   typedef typename Rhs::Scalar Scalar;
0068   typedef evaluator<Lhs> LhsEval;
0069   typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
0070   static void run(const Lhs& lhs, Rhs& other)
0071   {
0072     LhsEval lhsEval(lhs);
0073     for(Index col=0 ; col<other.cols() ; ++col)
0074     {
0075       for(Index i=lhs.rows()-1 ; i>=0 ; --i)
0076       {
0077         Scalar tmp = other.coeff(i,col);
0078         Scalar l_ii(0);
0079         LhsIterator it(lhsEval, i);
0080         while(it && it.index()<i)
0081           ++it;
0082         if(!(Mode & UnitDiag))
0083         {
0084           eigen_assert(it && it.index()==i);
0085           l_ii = it.value();
0086           ++it;
0087         }
0088         else if (it && it.index() == i)
0089           ++it;
0090         for(; it; ++it)
0091         {
0092           tmp -= it.value() * other.coeff(it.index(),col);
0093         }
0094
0095         if (Mode & UnitDiag)  other.coeffRef(i,col) = tmp;
0096         else                  other.coeffRef(i,col) = tmp/l_ii;
0097       }
0098     }
0099   }
0100 };
0101
0102 // forward substitution, col-major
0103 template<typename Lhs, typename Rhs, int Mode>
0104 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,ColMajor>
0105 {
0106   typedef typename Rhs::Scalar Scalar;
0107   typedef evaluator<Lhs> LhsEval;
0108   typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
0109   static void run(const Lhs& lhs, Rhs& other)
0110   {
0111     LhsEval lhsEval(lhs);
0112     for(Index col=0 ; col<other.cols() ; ++col)
0113     {
0114       for(Index i=0; i<lhs.cols(); ++i)
0115       {
0116         Scalar& tmp = other.coeffRef(i,col);
0117         if (tmp!=Scalar(0)) // optimization when other is actually sparse
0118         {
0119           LhsIterator it(lhsEval, i);
0120           while(it && it.index()<i)
0121             ++it;
0122           if(!(Mode & UnitDiag))
0123           {
0124             eigen_assert(it && it.index()==i);
0125             tmp /= it.value();
0126           }
0127           if (it && it.index()==i)
0128             ++it;
0129           for(; it; ++it)
0130             other.coeffRef(it.index(), col) -= tmp * it.value();
0131         }
0132       }
0133     }
0134   }
0135 };
0136
0137 // backward substitution, col-major
0138 template<typename Lhs, typename Rhs, int Mode>
0139 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,ColMajor>
0140 {
0141   typedef typename Rhs::Scalar Scalar;
0142   typedef evaluator<Lhs> LhsEval;
0143   typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
0144   static void run(const Lhs& lhs, Rhs& other)
0145   {
0146     LhsEval lhsEval(lhs);
0147     for(Index col=0 ; col<other.cols() ; ++col)
0148     {
0149       for(Index i=lhs.cols()-1; i>=0; --i)
0150       {
0151         Scalar& tmp = other.coeffRef(i,col);
0152         if (tmp!=Scalar(0)) // optimization when other is actually sparse
0153         {
0154           if(!(Mode & UnitDiag))
0155           {
0156             // TODO replace this by a binary search. make sure the binary search is safe for partially sorted elements
0157             LhsIterator it(lhsEval, i);
0158             while(it && it.index()!=i)
0159               ++it;
0160             eigen_assert(it && it.index()==i);
0161             other.coeffRef(i,col) /= it.value();
0162           }
0163           LhsIterator it(lhsEval, i);
0164           for(; it && it.index()<i; ++it)
0165             other.coeffRef(it.index(), col) -= tmp * it.value();
0166         }
0167       }
0168     }
0169   }
0170 };
0171
0172 } // end namespace internal
0173
0174 #ifndef EIGEN_PARSED_BY_DOXYGEN
0175
0176 template<typename ExpressionType,unsigned int Mode>
0177 template<typename OtherDerived>
0178 void TriangularViewImpl<ExpressionType,Mode,Sparse>::solveInPlace(MatrixBase<OtherDerived>& other) const
0179 {
0180   eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows());
0181   eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
0182
0183   enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
0184
0185   typedef typename internal::conditional<copy,
0186     typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
0187   OtherCopy otherCopy(other.derived());
0188
0189   internal::sparse_solve_triangular_selector<ExpressionType, typename internal::remove_reference<OtherCopy>::type, Mode>::run(derived().nestedExpression(), otherCopy);
0190
0191   if (copy)
0192     other = otherCopy;
0193 }
0194 #endif
0195
0196 // pure sparse path
0197
0198 namespace internal {
0199
0200 template<typename Lhs, typename Rhs, int Mode,
0201   int UpLo = (Mode & Lower)
0202            ? Lower
0203            : (Mode & Upper)
0204            ? Upper
0205            : -1,
0206   int StorageOrder = int(Lhs::Flags) & (RowMajorBit)>
0207 struct sparse_solve_triangular_sparse_selector;
0208
0209 // forward substitution, col-major
0210 template<typename Lhs, typename Rhs, int Mode, int UpLo>
0211 struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
0212 {
0213   typedef typename Rhs::Scalar Scalar;
0214   typedef typename promote_index_type<typename traits<Lhs>::StorageIndex,
0215                                       typename traits<Rhs>::StorageIndex>::type StorageIndex;
0216   static void run(const Lhs& lhs, Rhs& other)
0217   {
0218     const bool IsLower = (UpLo==Lower);
0219     AmbiVector<Scalar,StorageIndex> tempVector(other.rows()*2);
0220     tempVector.setBounds(0,other.rows());
0221
0222     Rhs res(other.rows(), other.cols());
0223     res.reserve(other.nonZeros());
0224
0225     for(Index col=0 ; col<other.cols() ; ++col)
0226     {
0227       // FIXME estimate number of non zeros
0228       tempVector.init(.99/*float(other.col(col).nonZeros())/float(other.rows())*/);
0229       tempVector.setZero();
0230       tempVector.restart();
0231       for (typename Rhs::InnerIterator rhsIt(other, col); rhsIt; ++rhsIt)
0232       {
0233         tempVector.coeffRef(rhsIt.index()) = rhsIt.value();
0234       }
0235
0236       for(Index i=IsLower?0:lhs.cols()-1;
0237           IsLower?i<lhs.cols():i>=0;
0238           i+=IsLower?1:-1)
0239       {
0240         tempVector.restart();
0241         Scalar& ci = tempVector.coeffRef(i);
0242         if (ci!=Scalar(0))
0243         {
0244           // find
0245           typename Lhs::InnerIterator it(lhs, i);
0246           if(!(Mode & UnitDiag))
0247           {
0248             if (IsLower)
0249             {
0250               eigen_assert(it.index()==i);
0251               ci /= it.value();
0252             }
0253             else
0254               ci /= lhs.coeff(i,i);
0255           }
0256           tempVector.restart();
0257           if (IsLower)
0258           {
0259             if (it.index()==i)
0260               ++it;
0261             for(; it; ++it)
0262               tempVector.coeffRef(it.index()) -= ci * it.value();
0263           }
0264           else
0265           {
0266             for(; it && it.index()<i; ++it)
0267               tempVector.coeffRef(it.index()) -= ci * it.value();
0268           }
0269         }
0270       }
0271
0272
0273       Index count = 0;
0274       // FIXME compute a reference value to filter zeros
0275       for (typename AmbiVector<Scalar,StorageIndex>::Iterator it(tempVector/*,1e-12*/); it; ++it)
0276       {
0277         ++ count;
0278 //         std::cerr << "fill " << it.index() << ", " << col << "\n";
0279 //         std::cout << it.value() << "  ";
0280         // FIXME use insertBack
0281         res.insert(it.index(), col) = it.value();
0282       }
0283 //       std::cout << "tempVector.nonZeros() == " << int(count) << " / " << (other.rows()) << "\n";
0284     }
0285     res.finalize();
0286     other = res.markAsRValue();
0287   }
0288 };
0289
0290 } // end namespace internal
0291
0292 #ifndef EIGEN_PARSED_BY_DOXYGEN
0293 template<typename ExpressionType,unsigned int Mode>
0294 template<typename OtherDerived>
0295 void TriangularViewImpl<ExpressionType,Mode,Sparse>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const
0296 {
0297   eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows());
0298   eigen_assert( (!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
0299
0300 //   enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
0301
0302 //   typedef typename internal::conditional<copy,
0303 //     typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
0304 //   OtherCopy otherCopy(other.derived());
0305
0306   internal::sparse_solve_triangular_sparse_selector<ExpressionType, OtherDerived, Mode>::run(derived().nestedExpression(), other.derived());
0307
0308 //   if (copy)
0309 //     other = otherCopy;
0310 }
0311 #endif
0312
0313 } // end namespace Eigen
0314
0315 #endif // EIGEN_SPARSETRIANGULARSOLVER_H