Eigen  3.3.3
TriangularSolver.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_SPARSETRIANGULARSOLVER_H
00011 #define EIGEN_SPARSETRIANGULARSOLVER_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal {
00016 
00017 template<typename Lhs, typename Rhs, int Mode,
00018   int UpLo = (Mode & Lower)
00019            ? Lower
00020            : (Mode & Upper)
00021            ? Upper
00022            : -1,
00023   int StorageOrder = int(traits<Lhs>::Flags) & RowMajorBit>
00024 struct sparse_solve_triangular_selector;
00025 
00026 // forward substitution, row-major
00027 template<typename Lhs, typename Rhs, int Mode>
00028 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,RowMajor>
00029 {
00030   typedef typename Rhs::Scalar Scalar;
00031   typedef evaluator<Lhs> LhsEval;
00032   typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
00033   static void run(const Lhs& lhs, Rhs& other)
00034   {
00035     LhsEval lhsEval(lhs);
00036     for(Index col=0 ; col<other.cols() ; ++col)
00037     {
00038       for(Index i=0; i<lhs.rows(); ++i)
00039       {
00040         Scalar tmp = other.coeff(i,col);
00041         Scalar lastVal(0);
00042         Index lastIndex = 0;
00043         for(LhsIterator it(lhsEval, i); it; ++it)
00044         {
00045           lastVal = it.value();
00046           lastIndex = it.index();
00047           if(lastIndex==i)
00048             break;
00049           tmp -= lastVal * other.coeff(lastIndex,col);
00050         }
00051         if (Mode & UnitDiag)
00052           other.coeffRef(i,col) = tmp;
00053         else
00054         {
00055           eigen_assert(lastIndex==i);
00056           other.coeffRef(i,col) = tmp/lastVal;
00057         }
00058       }
00059     }
00060   }
00061 };
00062 
00063 // backward substitution, row-major
00064 template<typename Lhs, typename Rhs, int Mode>
00065 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,RowMajor>
00066 {
00067   typedef typename Rhs::Scalar Scalar;
00068   typedef evaluator<Lhs> LhsEval;
00069   typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
00070   static void run(const Lhs& lhs, Rhs& other)
00071   {
00072     LhsEval lhsEval(lhs);
00073     for(Index col=0 ; col<other.cols() ; ++col)
00074     {
00075       for(Index i=lhs.rows()-1 ; i>=0 ; --i)
00076       {
00077         Scalar tmp = other.coeff(i,col);
00078         Scalar l_ii(0);
00079         LhsIterator it(lhsEval, i);
00080         while(it && it.index()<i)
00081           ++it;
00082         if(!(Mode & UnitDiag))
00083         {
00084           eigen_assert(it && it.index()==i);
00085           l_ii = it.value();
00086           ++it;
00087         }
00088         else if (it && it.index() == i)
00089           ++it;
00090         for(; it; ++it)
00091         {
00092           tmp -= it.value() * other.coeff(it.index(),col);
00093         }
00094 
00095         if (Mode & UnitDiag)  other.coeffRef(i,col) = tmp;
00096         else                  other.coeffRef(i,col) = tmp/l_ii;
00097       }
00098     }
00099   }
00100 };
00101 
00102 // forward substitution, col-major
00103 template<typename Lhs, typename Rhs, int Mode>
00104 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,ColMajor>
00105 {
00106   typedef typename Rhs::Scalar Scalar;
00107   typedef evaluator<Lhs> LhsEval;
00108   typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
00109   static void run(const Lhs& lhs, Rhs& other)
00110   {
00111     LhsEval lhsEval(lhs);
00112     for(Index col=0 ; col<other.cols() ; ++col)
00113     {
00114       for(Index i=0; i<lhs.cols(); ++i)
00115       {
00116         Scalar& tmp = other.coeffRef(i,col);
00117         if (tmp!=Scalar(0)) // optimization when other is actually sparse
00118         {
00119           LhsIterator it(lhsEval, i);
00120           while(it && it.index()<i)
00121             ++it;
00122           if(!(Mode & UnitDiag))
00123           {
00124             eigen_assert(it && it.index()==i);
00125             tmp /= it.value();
00126           }
00127           if (it && it.index()==i)
00128             ++it;
00129           for(; it; ++it)
00130             other.coeffRef(it.index(), col) -= tmp * it.value();
00131         }
00132       }
00133     }
00134   }
00135 };
00136 
00137 // backward substitution, col-major
00138 template<typename Lhs, typename Rhs, int Mode>
00139 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,ColMajor>
00140 {
00141   typedef typename Rhs::Scalar Scalar;
00142   typedef evaluator<Lhs> LhsEval;
00143   typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
00144   static void run(const Lhs& lhs, Rhs& other)
00145   {
00146     LhsEval lhsEval(lhs);
00147     for(Index col=0 ; col<other.cols() ; ++col)
00148     {
00149       for(Index i=lhs.cols()-1; i>=0; --i)
00150       {
00151         Scalar& tmp = other.coeffRef(i,col);
00152         if (tmp!=Scalar(0)) // optimization when other is actually sparse
00153         {
00154           if(!(Mode & UnitDiag))
00155           {
00156             // TODO replace this by a binary search. make sure the binary search is safe for partially sorted elements
00157             LhsIterator it(lhsEval, i);
00158             while(it && it.index()!=i)
00159               ++it;
00160             eigen_assert(it && it.index()==i);
00161             other.coeffRef(i,col) /= it.value();
00162           }
00163           LhsIterator it(lhsEval, i);
00164           for(; it && it.index()<i; ++it)
00165             other.coeffRef(it.index(), col) -= tmp * it.value();
00166         }
00167       }
00168     }
00169   }
00170 };
00171 
00172 } // end namespace internal
00173 
00174 #ifndef EIGEN_PARSED_BY_DOXYGEN
00175 
00176 template<typename ExpressionType,unsigned int Mode>
00177 template<typename OtherDerived>
00178 void TriangularViewImpl<ExpressionType,Mode,Sparse>::solveInPlace(MatrixBase<OtherDerived>& other) const
00179 {
00180   eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows());
00181   eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
00182 
00183   enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
00184 
00185   typedef typename internal::conditional<copy,
00186     typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
00187   OtherCopy otherCopy(other.derived());
00188 
00189   internal::sparse_solve_triangular_selector<ExpressionType, typename internal::remove_reference<OtherCopy>::type, Mode>::run(derived().nestedExpression(), otherCopy);
00190 
00191   if (copy)
00192     other = otherCopy;
00193 }
00194 #endif
00195 
00196 // pure sparse path
00197 
00198 namespace internal {
00199 
00200 template<typename Lhs, typename Rhs, int Mode,
00201   int UpLo = (Mode & Lower)
00202            ? Lower
00203            : (Mode & Upper)
00204            ? Upper
00205            : -1,
00206   int StorageOrder = int(Lhs::Flags) & (RowMajorBit)>
00207 struct sparse_solve_triangular_sparse_selector;
00208 
00209 // forward substitution, col-major
00210 template<typename Lhs, typename Rhs, int Mode, int UpLo>
00211 struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
00212 {
00213   typedef typename Rhs::Scalar Scalar;
00214   typedef typename promote_index_type<typename traits<Lhs>::StorageIndex,
00215                                       typename traits<Rhs>::StorageIndex>::type StorageIndex;
00216   static void run(const Lhs& lhs, Rhs& other)
00217   {
00218     const bool IsLower = (UpLo==Lower);
00219     AmbiVector<Scalar,StorageIndex> tempVector(other.rows()*2);
00220     tempVector.setBounds(0,other.rows());
00221 
00222     Rhs res(other.rows(), other.cols());
00223     res.reserve(other.nonZeros());
00224 
00225     for(Index col=0 ; col<other.cols() ; ++col)
00226     {
00227       // FIXME estimate number of non zeros
00228       tempVector.init(.99/*float(other.col(col).nonZeros())/float(other.rows())*/);
00229       tempVector.setZero();
00230       tempVector.restart();
00231       for (typename Rhs::InnerIterator rhsIt(other, col); rhsIt; ++rhsIt)
00232       {
00233         tempVector.coeffRef(rhsIt.index()) = rhsIt.value();
00234       }
00235 
00236       for(Index i=IsLower?0:lhs.cols()-1;
00237           IsLower?i<lhs.cols():i>=0;
00238           i+=IsLower?1:-1)
00239       {
00240         tempVector.restart();
00241         Scalar& ci = tempVector.coeffRef(i);
00242         if (ci!=Scalar(0))
00243         {
00244           // find
00245           typename Lhs::InnerIterator it(lhs, i);
00246           if(!(Mode & UnitDiag))
00247           {
00248             if (IsLower)
00249             {
00250               eigen_assert(it.index()==i);
00251               ci /= it.value();
00252             }
00253             else
00254               ci /= lhs.coeff(i,i);
00255           }
00256           tempVector.restart();
00257           if (IsLower)
00258           {
00259             if (it.index()==i)
00260               ++it;
00261             for(; it; ++it)
00262               tempVector.coeffRef(it.index()) -= ci * it.value();
00263           }
00264           else
00265           {
00266             for(; it && it.index()<i; ++it)
00267               tempVector.coeffRef(it.index()) -= ci * it.value();
00268           }
00269         }
00270       }
00271 
00272 
00273       Index count = 0;
00274       // FIXME compute a reference value to filter zeros
00275       for (typename AmbiVector<Scalar,StorageIndex>::Iterator it(tempVector/*,1e-12*/); it; ++it)
00276       {
00277         ++ count;
00278 //         std::cerr << "fill " << it.index() << ", " << col << "\n";
00279 //         std::cout << it.value() << "  ";
00280         // FIXME use insertBack
00281         res.insert(it.index(), col) = it.value();
00282       }
00283 //       std::cout << "tempVector.nonZeros() == " << int(count) << " / " << (other.rows()) << "\n";
00284     }
00285     res.finalize();
00286     other = res.markAsRValue();
00287   }
00288 };
00289 
00290 } // end namespace internal
00291 
00292 #ifndef EIGEN_PARSED_BY_DOXYGEN
00293 template<typename ExpressionType,unsigned int Mode>
00294 template<typename OtherDerived>
00295 void TriangularViewImpl<ExpressionType,Mode,Sparse>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const
00296 {
00297   eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows());
00298   eigen_assert( (!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
00299 
00300 //   enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
00301 
00302 //   typedef typename internal::conditional<copy,
00303 //     typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
00304 //   OtherCopy otherCopy(other.derived());
00305 
00306   internal::sparse_solve_triangular_sparse_selector<ExpressionType, OtherDerived, Mode>::run(derived().nestedExpression(), other.derived());
00307 
00308 //   if (copy)
00309 //     other = otherCopy;
00310 }
00311 #endif
00312 
00313 } // end namespace Eigen
00314 
00315 #endif // EIGEN_SPARSETRIANGULARSOLVER_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends