$treeview $search $mathjax
Eigen-unsupported
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) 2012 Giacomo Po <gpo@ucla.edu> 00005 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 00012 #ifndef EIGEN_MINRES_H_ 00013 #define EIGEN_MINRES_H_ 00014 00015 00016 namespace Eigen { 00017 00018 namespace internal { 00019 00029 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner> 00030 EIGEN_DONT_INLINE 00031 void minres(const MatrixType& mat, const Rhs& rhs, Dest& x, 00032 const Preconditioner& precond, int& iters, 00033 typename Dest::RealScalar& tol_error) 00034 { 00035 using std::sqrt; 00036 typedef typename Dest::RealScalar RealScalar; 00037 typedef typename Dest::Scalar Scalar; 00038 typedef Matrix<Scalar,Dynamic,1> VectorType; 00039 00040 // Check for zero rhs 00041 const RealScalar rhsNorm2(rhs.squaredNorm()); 00042 if(rhsNorm2 == 0) 00043 { 00044 x.setZero(); 00045 iters = 0; 00046 tol_error = 0; 00047 return; 00048 } 00049 00050 // initialize 00051 const int maxIters(iters); // initialize maxIters to iters 00052 const int N(mat.cols()); // the size of the matrix 00053 const RealScalar threshold2(tol_error*tol_error*rhsNorm2); // convergence threshold (compared to residualNorm2) 00054 00055 // Initialize preconditioned Lanczos 00056 VectorType v_old(N); // will be initialized inside loop 00057 VectorType v( VectorType::Zero(N) ); //initialize v 00058 VectorType v_new(rhs-mat*x); //initialize v_new 00059 RealScalar residualNorm2(v_new.squaredNorm()); 00060 VectorType w(N); // will be initialized inside loop 00061 VectorType w_new(precond.solve(v_new)); // initialize w_new 00062 // RealScalar beta; // will be initialized inside loop 00063 RealScalar beta_new2(v_new.dot(w_new)); 00064 eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE"); 00065 RealScalar beta_new(sqrt(beta_new2)); 00066 const RealScalar beta_one(beta_new); 00067 v_new /= beta_new; 00068 w_new /= beta_new; 00069 // Initialize other variables 00070 RealScalar c(1.0); // the cosine of the Givens rotation 00071 RealScalar c_old(1.0); 00072 RealScalar s(0.0); // the sine of the Givens rotation 00073 RealScalar s_old(0.0); // the sine of the Givens rotation 00074 VectorType p_oold(N); // will be initialized in loop 00075 VectorType p_old(VectorType::Zero(N)); // initialize p_old=0 00076 VectorType p(p_old); // initialize p=0 00077 RealScalar eta(1.0); 00078 00079 iters = 0; // reset iters 00080 while ( iters < maxIters ) 00081 { 00082 // Preconditioned Lanczos 00083 /* Note that there are 4 variants on the Lanczos algorithm. These are 00084 * described in Paige, C. C. (1972). Computational variants of 00085 * the Lanczos method for the eigenproblem. IMA Journal of Applied 00086 * Mathematics, 10(3), 373–381. The current implementation corresponds 00087 * to the case A(2,7) in the paper. It also corresponds to 00088 * algorithm 6.14 in Y. Saad, Iterative Methods for Sparse Linear 00089 * Systems, 2003 p.173. For the preconditioned version see 00090 * A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM (1987). 00091 */ 00092 const RealScalar beta(beta_new); 00093 v_old = v; // update: at first time step, this makes v_old = 0 so value of beta doesn't matter 00094 // const VectorType v_old(v); // NOT SURE IF CREATING v_old EVERY ITERATION IS EFFICIENT 00095 v = v_new; // update 00096 w = w_new; // update 00097 // const VectorType w(w_new); // NOT SURE IF CREATING w EVERY ITERATION IS EFFICIENT 00098 v_new.noalias() = mat*w - beta*v_old; // compute v_new 00099 const RealScalar alpha = v_new.dot(w); 00100 v_new -= alpha*v; // overwrite v_new 00101 w_new = precond.solve(v_new); // overwrite w_new 00102 beta_new2 = v_new.dot(w_new); // compute beta_new 00103 eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE"); 00104 beta_new = sqrt(beta_new2); // compute beta_new 00105 v_new /= beta_new; // overwrite v_new for next iteration 00106 w_new /= beta_new; // overwrite w_new for next iteration 00107 00108 // Givens rotation 00109 const RealScalar r2 =s*alpha+c*c_old*beta; // s, s_old, c and c_old are still from previous iteration 00110 const RealScalar r3 =s_old*beta; // s, s_old, c and c_old are still from previous iteration 00111 const RealScalar r1_hat=c*alpha-c_old*s*beta; 00112 const RealScalar r1 =sqrt( std::pow(r1_hat,2) + std::pow(beta_new,2) ); 00113 c_old = c; // store for next iteration 00114 s_old = s; // store for next iteration 00115 c=r1_hat/r1; // new cosine 00116 s=beta_new/r1; // new sine 00117 00118 // Update solution 00119 p_oold = p_old; 00120 // const VectorType p_oold(p_old); // NOT SURE IF CREATING p_oold EVERY ITERATION IS EFFICIENT 00121 p_old = p; 00122 p.noalias()=(w-r2*p_old-r3*p_oold) /r1; // IS NOALIAS REQUIRED? 00123 x += beta_one*c*eta*p; 00124 00125 /* Update the squared residual. Note that this is the estimated residual. 00126 The real residual |Ax-b|^2 may be slightly larger */ 00127 residualNorm2 *= s*s; 00128 00129 if ( residualNorm2 < threshold2) 00130 { 00131 break; 00132 } 00133 00134 eta=-s*eta; // update eta 00135 iters++; // increment iteration number (for output purposes) 00136 } 00137 00138 /* Compute error. Note that this is the estimated error. The real 00139 error |Ax-b|/|b| may be slightly larger */ 00140 tol_error = std::sqrt(residualNorm2 / rhsNorm2); 00141 } 00142 00143 } 00144 00145 template< typename _MatrixType, int _UpLo=Lower, 00146 typename _Preconditioner = IdentityPreconditioner> 00147 // typename _Preconditioner = IdentityPreconditioner<typename _MatrixType::Scalar> > // preconditioner must be positive definite 00148 class MINRES; 00149 00150 namespace internal { 00151 00152 template< typename _MatrixType, int _UpLo, typename _Preconditioner> 00153 struct traits<MINRES<_MatrixType,_UpLo,_Preconditioner> > 00154 { 00155 typedef _MatrixType MatrixType; 00156 typedef _Preconditioner Preconditioner; 00157 }; 00158 00159 } 00160 00197 template< typename _MatrixType, int _UpLo, typename _Preconditioner> 00198 class MINRES : public IterativeSolverBase<MINRES<_MatrixType,_UpLo,_Preconditioner> > 00199 { 00200 00201 typedef IterativeSolverBase<MINRES> Base; 00202 using Base::mp_matrix; 00203 using Base::m_error; 00204 using Base::m_iterations; 00205 using Base::m_info; 00206 using Base::m_isInitialized; 00207 public: 00208 typedef _MatrixType MatrixType; 00209 typedef typename MatrixType::Scalar Scalar; 00210 typedef typename MatrixType::Index Index; 00211 typedef typename MatrixType::RealScalar RealScalar; 00212 typedef _Preconditioner Preconditioner; 00213 00214 enum {UpLo = _UpLo}; 00215 00216 public: 00217 00219 MINRES() : Base() {} 00220 00231 MINRES(const MatrixType& A) : Base(A) {} 00232 00234 ~MINRES(){} 00235 00241 template<typename Rhs,typename Guess> 00242 inline const internal::solve_retval_with_guess<MINRES, Rhs, Guess> 00243 solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const 00244 { 00245 eigen_assert(m_isInitialized && "MINRES is not initialized."); 00246 eigen_assert(Base::rows()==b.rows() 00247 && "MINRES::solve(): invalid number of rows of the right hand side matrix b"); 00248 return internal::solve_retval_with_guess 00249 <MINRES, Rhs, Guess>(*this, b.derived(), x0); 00250 } 00251 00253 template<typename Rhs,typename Dest> 00254 void _solveWithGuess(const Rhs& b, Dest& x) const 00255 { 00256 typedef typename internal::conditional<UpLo==(Lower|Upper), 00257 const MatrixType&, 00258 SparseSelfAdjointView<const MatrixType, UpLo> 00259 >::type MatrixWrapperType; 00260 00261 m_iterations = Base::maxIterations(); 00262 m_error = Base::m_tolerance; 00263 00264 for(int j=0; j<b.cols(); ++j) 00265 { 00266 m_iterations = Base::maxIterations(); 00267 m_error = Base::m_tolerance; 00268 00269 typename Dest::ColXpr xj(x,j); 00270 internal::minres(MatrixWrapperType(*mp_matrix), b.col(j), xj, 00271 Base::m_preconditioner, m_iterations, m_error); 00272 } 00273 00274 m_isInitialized = true; 00275 m_info = m_error <= Base::m_tolerance ? Success : NoConvergence; 00276 } 00277 00279 template<typename Rhs,typename Dest> 00280 void _solve(const Rhs& b, Dest& x) const 00281 { 00282 x.setZero(); 00283 _solveWithGuess(b,x); 00284 } 00285 00286 protected: 00287 00288 }; 00289 00290 namespace internal { 00291 00292 template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs> 00293 struct solve_retval<MINRES<_MatrixType,_UpLo,_Preconditioner>, Rhs> 00294 : solve_retval_base<MINRES<_MatrixType,_UpLo,_Preconditioner>, Rhs> 00295 { 00296 typedef MINRES<_MatrixType,_UpLo,_Preconditioner> Dec; 00297 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 00298 00299 template<typename Dest> void evalTo(Dest& dst) const 00300 { 00301 dec()._solve(rhs(),dst); 00302 } 00303 }; 00304 00305 } // end namespace internal 00306 00307 } // end namespace Eigen 00308 00309 #endif // EIGEN_MINRES_H 00310