![]() |
Eigen
3.3.3
|
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