![]() |
Eigen-unsupported
3.3.3
|
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-2014 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, Index& 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 Index maxIters(iters); // initialize maxIters to iters 00052 const Index 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 class MINRES; 00148 00149 namespace internal { 00150 00151 template< typename _MatrixType, int _UpLo, typename _Preconditioner> 00152 struct traits<MINRES<_MatrixType,_UpLo,_Preconditioner> > 00153 { 00154 typedef _MatrixType MatrixType; 00155 typedef _Preconditioner Preconditioner; 00156 }; 00157 00158 } 00159 00198 template< typename _MatrixType, int _UpLo, typename _Preconditioner> 00199 class MINRES : public IterativeSolverBase<MINRES<_MatrixType,_UpLo,_Preconditioner> > 00200 { 00201 00202 typedef IterativeSolverBase<MINRES> Base; 00203 using Base::matrix; 00204 using Base::m_error; 00205 using Base::m_iterations; 00206 using Base::m_info; 00207 using Base::m_isInitialized; 00208 public: 00209 using Base::_solve_impl; 00210 typedef _MatrixType MatrixType; 00211 typedef typename MatrixType::Scalar Scalar; 00212 typedef typename MatrixType::RealScalar RealScalar; 00213 typedef _Preconditioner Preconditioner; 00214 00215 enum {UpLo = _UpLo}; 00216 00217 public: 00218 00220 MINRES() : Base() {} 00221 00232 template<typename MatrixDerived> 00233 explicit MINRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {} 00234 00236 ~MINRES(){} 00237 00239 template<typename Rhs,typename Dest> 00240 void _solve_with_guess_impl(const Rhs& b, Dest& x) const 00241 { 00242 typedef typename Base::MatrixWrapper MatrixWrapper; 00243 typedef typename Base::ActualMatrixType ActualMatrixType; 00244 enum { 00245 TransposeInput = (!MatrixWrapper::MatrixFree) 00246 && (UpLo==(Lower|Upper)) 00247 && (!MatrixType::IsRowMajor) 00248 && (!NumTraits<Scalar>::IsComplex) 00249 }; 00250 typedef typename internal::conditional<TransposeInput,Transpose<const ActualMatrixType>, ActualMatrixType const&>::type RowMajorWrapper; 00251 EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(MatrixWrapper::MatrixFree,UpLo==(Lower|Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY); 00252 typedef typename internal::conditional<UpLo==(Lower|Upper), 00253 RowMajorWrapper, 00254 typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type 00255 >::type SelfAdjointWrapper; 00256 00257 m_iterations = Base::maxIterations(); 00258 m_error = Base::m_tolerance; 00259 RowMajorWrapper row_mat(matrix()); 00260 for(int j=0; j<b.cols(); ++j) 00261 { 00262 m_iterations = Base::maxIterations(); 00263 m_error = Base::m_tolerance; 00264 00265 typename Dest::ColXpr xj(x,j); 00266 internal::minres(SelfAdjointWrapper(row_mat), b.col(j), xj, 00267 Base::m_preconditioner, m_iterations, m_error); 00268 } 00269 00270 m_isInitialized = true; 00271 m_info = m_error <= Base::m_tolerance ? Success : NoConvergence; 00272 } 00273 00275 template<typename Rhs,typename Dest> 00276 void _solve_impl(const Rhs& b, MatrixBase<Dest> &x) const 00277 { 00278 x.setZero(); 00279 _solve_with_guess_impl(b,x.derived()); 00280 } 00281 00282 protected: 00283 00284 }; 00285 00286 } // end namespace Eigen 00287 00288 #endif // EIGEN_MINRES_H 00289