$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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