TriangularSolver.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 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   static void run(const Lhs& lhs, Rhs& other)
00032   {
00033     for(int col=0 ; col<other.cols() ; ++col)
00034     {
00035       for(int i=0; i<lhs.rows(); ++i)
00036       {
00037         Scalar tmp = other.coeff(i,col);
00038         Scalar lastVal(0);
00039         int lastIndex = 0;
00040         for(typename Lhs::InnerIterator it(lhs, i); it; ++it)
00041         {
00042           lastVal = it.value();
00043           lastIndex = it.index();
00044           if(lastIndex==i)
00045             break;
00046           tmp -= lastVal * other.coeff(lastIndex,col);
00047         }
00048         if (Mode & UnitDiag)
00049           other.coeffRef(i,col) = tmp;
00050         else
00051         {
00052           eigen_assert(lastIndex==i);
00053           other.coeffRef(i,col) = tmp/lastVal;
00054         }
00055       }
00056     }
00057   }
00058 };
00059 
00060 // backward substitution, row-major
00061 template<typename Lhs, typename Rhs, int Mode>
00062 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,RowMajor>
00063 {
00064   typedef typename Rhs::Scalar Scalar;
00065   static void run(const Lhs& lhs, Rhs& other)
00066   {
00067     for(int col=0 ; col<other.cols() ; ++col)
00068     {
00069       for(int i=lhs.rows()-1 ; i>=0 ; --i)
00070       {
00071         Scalar tmp = other.coeff(i,col);
00072         Scalar l_ii = 0;
00073         typename Lhs::InnerIterator it(lhs, i);
00074         while(it && it.index()<i)
00075           ++it;
00076         if(!(Mode & UnitDiag))
00077         {
00078           eigen_assert(it && it.index()==i);
00079           l_ii = it.value();
00080           ++it;
00081         }
00082         else if (it && it.index() == i)
00083           ++it;
00084         for(; it; ++it)
00085         {
00086           tmp -= it.value() * other.coeff(it.index(),col);
00087         }
00088 
00089         if (Mode & UnitDiag)
00090           other.coeffRef(i,col) = tmp;
00091         else
00092           other.coeffRef(i,col) = tmp/l_ii;
00093       }
00094     }
00095   }
00096 };
00097 
00098 // forward substitution, col-major
00099 template<typename Lhs, typename Rhs, int Mode>
00100 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,ColMajor>
00101 {
00102   typedef typename Rhs::Scalar Scalar;
00103   static void run(const Lhs& lhs, Rhs& other)
00104   {
00105     for(int col=0 ; col<other.cols() ; ++col)
00106     {
00107       for(int i=0; i<lhs.cols(); ++i)
00108       {
00109         Scalar& tmp = other.coeffRef(i,col);
00110         if (tmp!=Scalar(0)) // optimization when other is actually sparse
00111         {
00112           typename Lhs::InnerIterator it(lhs, i);
00113           while(it && it.index()<i)
00114             ++it;
00115           if(!(Mode & UnitDiag))
00116           {
00117             eigen_assert(it && it.index()==i);
00118             tmp /= it.value();
00119           }
00120           if (it && it.index()==i)
00121             ++it;
00122           for(; it; ++it)
00123             other.coeffRef(it.index(), col) -= tmp * it.value();
00124         }
00125       }
00126     }
00127   }
00128 };
00129 
00130 // backward substitution, col-major
00131 template<typename Lhs, typename Rhs, int Mode>
00132 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,ColMajor>
00133 {
00134   typedef typename Rhs::Scalar Scalar;
00135   static void run(const Lhs& lhs, Rhs& other)
00136   {
00137     for(int col=0 ; col<other.cols() ; ++col)
00138     {
00139       for(int i=lhs.cols()-1; i>=0; --i)
00140       {
00141         Scalar& tmp = other.coeffRef(i,col);
00142         if (tmp!=Scalar(0)) // optimization when other is actually sparse
00143         {
00144           if(!(Mode & UnitDiag))
00145           {
00146             // TODO replace this by a binary search. make sure the binary search is safe for partially sorted elements
00147             typename Lhs::ReverseInnerIterator it(lhs, i);
00148             while(it && it.index()!=i)
00149               --it;
00150             eigen_assert(it && it.index()==i);
00151             other.coeffRef(i,col) /= it.value();
00152           }
00153           typename Lhs::InnerIterator it(lhs, i);
00154           for(; it && it.index()<i; ++it)
00155             other.coeffRef(it.index(), col) -= tmp * it.value();
00156         }
00157       }
00158     }
00159   }
00160 };
00161 
00162 } // end namespace internal
00163 
00164 template<typename ExpressionType,int Mode>
00165 template<typename OtherDerived>
00166 void SparseTriangularView<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDerived>& other) const
00167 {
00168   eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows());
00169   eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
00170 
00171   enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
00172 
00173   typedef typename internal::conditional<copy,
00174     typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
00175   OtherCopy otherCopy(other.derived());
00176 
00177   internal::sparse_solve_triangular_selector<ExpressionType, typename internal::remove_reference<OtherCopy>::type, Mode>::run(m_matrix, otherCopy);
00178 
00179   if (copy)
00180     other = otherCopy;
00181 }
00182 
00183 template<typename ExpressionType,int Mode>
00184 template<typename OtherDerived>
00185 typename internal::plain_matrix_type_column_major<OtherDerived>::type
00186 SparseTriangularView<ExpressionType,Mode>::solve(const MatrixBase<OtherDerived>& other) const
00187 {
00188   typename internal::plain_matrix_type_column_major<OtherDerived>::type res(other);
00189   solveInPlace(res);
00190   return res;
00191 }
00192 
00193 // pure sparse path
00194 
00195 namespace internal {
00196 
00197 template<typename Lhs, typename Rhs, int Mode,
00198   int UpLo = (Mode & Lower)
00199            ? Lower
00200            : (Mode & Upper)
00201            ? Upper
00202            : -1,
00203   int StorageOrder = int(Lhs::Flags) & (RowMajorBit)>
00204 struct sparse_solve_triangular_sparse_selector;
00205 
00206 // forward substitution, col-major
00207 template<typename Lhs, typename Rhs, int Mode, int UpLo>
00208 struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
00209 {
00210   typedef typename Rhs::Scalar Scalar;
00211   typedef typename promote_index_type<typename traits<Lhs>::Index,
00212                                          typename traits<Rhs>::Index>::type Index;
00213   static void run(const Lhs& lhs, Rhs& other)
00214   {
00215     const bool IsLower = (UpLo==Lower);
00216     AmbiVector<Scalar,Index> tempVector(other.rows()*2);
00217     tempVector.setBounds(0,other.rows());
00218 
00219     Rhs res(other.rows(), other.cols());
00220     res.reserve(other.nonZeros());
00221 
00222     for(int col=0 ; col<other.cols() ; ++col)
00223     {
00224       // FIXME estimate number of non zeros
00225       tempVector.init(.99/*float(other.col(col).nonZeros())/float(other.rows())*/);
00226       tempVector.setZero();
00227       tempVector.restart();
00228       for (typename Rhs::InnerIterator rhsIt(other, col); rhsIt; ++rhsIt)
00229       {
00230         tempVector.coeffRef(rhsIt.index()) = rhsIt.value();
00231       }
00232 
00233       for(int i=IsLower?0:lhs.cols()-1;
00234           IsLower?i<lhs.cols():i>=0;
00235           i+=IsLower?1:-1)
00236       {
00237         tempVector.restart();
00238         Scalar& ci = tempVector.coeffRef(i);
00239         if (ci!=Scalar(0))
00240         {
00241           // find
00242           typename Lhs::InnerIterator it(lhs, i);
00243           if(!(Mode & UnitDiag))
00244           {
00245             if (IsLower)
00246             {
00247               eigen_assert(it.index()==i);
00248               ci /= it.value();
00249             }
00250             else
00251               ci /= lhs.coeff(i,i);
00252           }
00253           tempVector.restart();
00254           if (IsLower)
00255           {
00256             if (it.index()==i)
00257               ++it;
00258             for(; it; ++it)
00259               tempVector.coeffRef(it.index()) -= ci * it.value();
00260           }
00261           else
00262           {
00263             for(; it && it.index()<i; ++it)
00264               tempVector.coeffRef(it.index()) -= ci * it.value();
00265           }
00266         }
00267       }
00268 
00269 
00270       int count = 0;
00271       // FIXME compute a reference value to filter zeros
00272       for (typename AmbiVector<Scalar,Index>::Iterator it(tempVector/*,1e-12*/); it; ++it)
00273       {
00274         ++ count;
00275 //         std::cerr << "fill " << it.index() << ", " << col << "\n";
00276 //         std::cout << it.value() << "  ";
00277         // FIXME use insertBack
00278         res.insert(it.index(), col) = it.value();
00279       }
00280 //       std::cout << "tempVector.nonZeros() == " << int(count) << " / " << (other.rows()) << "\n";
00281     }
00282     res.finalize();
00283     other = res.markAsRValue();
00284   }
00285 };
00286 
00287 } // end namespace internal
00288 
00289 template<typename ExpressionType,int Mode>
00290 template<typename OtherDerived>
00291 void SparseTriangularView<ExpressionType,Mode>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const
00292 {
00293   eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows());
00294   eigen_assert( (!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
00295 
00296 //   enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
00297 
00298 //   typedef typename internal::conditional<copy,
00299 //     typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
00300 //   OtherCopy otherCopy(other.derived());
00301 
00302   internal::sparse_solve_triangular_sparse_selector<ExpressionType, OtherDerived, Mode>::run(m_matrix, other.derived());
00303 
00304 //   if (copy)
00305 //     other = otherCopy;
00306 }
00307 
00308 #ifdef EIGEN2_SUPPORT
00309 
00310 // deprecated stuff:
00311 
00313 template<typename Derived>
00314 template<typename OtherDerived>
00315 void SparseMatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& other) const
00316 {
00317   this->template triangular<Flags&(Upper|Lower)>().solveInPlace(other);
00318 }
00319 
00321 template<typename Derived>
00322 template<typename OtherDerived>
00323 typename internal::plain_matrix_type_column_major<OtherDerived>::type
00324 SparseMatrixBase<Derived>::solveTriangular(const MatrixBase<OtherDerived>& other) const
00325 {
00326   typename internal::plain_matrix_type_column_major<OtherDerived>::type res(other);
00327   derived().solveTriangularInPlace(res);
00328   return res;
00329 }
00330 #endif // EIGEN2_SUPPORT
00331 
00332 } // end namespace Eigen
00333 
00334 #endif // EIGEN_SPARSETRIANGULARSOLVER_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 12:01:17