![]() |
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) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // Copyright (C) 2012, 2014 Kolja Brix <brix@igpm.rwth-aaachen.de> 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 #ifndef EIGEN_GMRES_H 00012 #define EIGEN_GMRES_H 00013 00014 namespace Eigen { 00015 00016 namespace internal { 00017 00055 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner> 00056 bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Preconditioner & precond, 00057 Index &iters, const Index &restart, typename Dest::RealScalar & tol_error) { 00058 00059 using std::sqrt; 00060 using std::abs; 00061 00062 typedef typename Dest::RealScalar RealScalar; 00063 typedef typename Dest::Scalar Scalar; 00064 typedef Matrix < Scalar, Dynamic, 1 > VectorType; 00065 typedef Matrix < Scalar, Dynamic, Dynamic, ColMajor> FMatrixType; 00066 00067 RealScalar tol = tol_error; 00068 const Index maxIters = iters; 00069 iters = 0; 00070 00071 const Index m = mat.rows(); 00072 00073 // residual and preconditioned residual 00074 VectorType p0 = rhs - mat*x; 00075 VectorType r0 = precond.solve(p0); 00076 00077 const RealScalar r0Norm = r0.norm(); 00078 00079 // is initial guess already good enough? 00080 if(r0Norm == 0) 00081 { 00082 tol_error = 0; 00083 return true; 00084 } 00085 00086 // storage for Hessenberg matrix and Householder data 00087 FMatrixType H = FMatrixType::Zero(m, restart + 1); 00088 VectorType w = VectorType::Zero(restart + 1); 00089 VectorType tau = VectorType::Zero(restart + 1); 00090 00091 // storage for Jacobi rotations 00092 std::vector < JacobiRotation < Scalar > > G(restart); 00093 00094 // storage for temporaries 00095 VectorType t(m), v(m), workspace(m), x_new(m); 00096 00097 // generate first Householder vector 00098 Ref<VectorType> H0_tail = H.col(0).tail(m - 1); 00099 RealScalar beta; 00100 r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta); 00101 w(0) = Scalar(beta); 00102 00103 for (Index k = 1; k <= restart; ++k) 00104 { 00105 ++iters; 00106 00107 v = VectorType::Unit(m, k - 1); 00108 00109 // apply Householder reflections H_{1} ... H_{k-1} to v 00110 // TODO: use a HouseholderSequence 00111 for (Index i = k - 1; i >= 0; --i) { 00112 v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); 00113 } 00114 00115 // apply matrix M to v: v = mat * v; 00116 t.noalias() = mat * v; 00117 v = precond.solve(t); 00118 00119 // apply Householder reflections H_{k-1} ... H_{1} to v 00120 // TODO: use a HouseholderSequence 00121 for (Index i = 0; i < k; ++i) { 00122 v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); 00123 } 00124 00125 if (v.tail(m - k).norm() != 0.0) 00126 { 00127 if (k <= restart) 00128 { 00129 // generate new Householder vector 00130 Ref<VectorType> Hk_tail = H.col(k).tail(m - k - 1); 00131 v.tail(m - k).makeHouseholder(Hk_tail, tau.coeffRef(k), beta); 00132 00133 // apply Householder reflection H_{k} to v 00134 v.tail(m - k).applyHouseholderOnTheLeft(Hk_tail, tau.coeffRef(k), workspace.data()); 00135 } 00136 } 00137 00138 if (k > 1) 00139 { 00140 for (Index i = 0; i < k - 1; ++i) 00141 { 00142 // apply old Givens rotations to v 00143 v.applyOnTheLeft(i, i + 1, G[i].adjoint()); 00144 } 00145 } 00146 00147 if (k<m && v(k) != (Scalar) 0) 00148 { 00149 // determine next Givens rotation 00150 G[k - 1].makeGivens(v(k - 1), v(k)); 00151 00152 // apply Givens rotation to v and w 00153 v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint()); 00154 w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint()); 00155 } 00156 00157 // insert coefficients into upper matrix triangle 00158 H.col(k-1).head(k) = v.head(k); 00159 00160 tol_error = abs(w(k)) / r0Norm; 00161 bool stop = (k==m || tol_error < tol || iters == maxIters); 00162 00163 if (stop || k == restart) 00164 { 00165 // solve upper triangular system 00166 Ref<VectorType> y = w.head(k); 00167 H.topLeftCorner(k, k).template triangularView <Upper>().solveInPlace(y); 00168 00169 // use Horner-like scheme to calculate solution vector 00170 x_new.setZero(); 00171 for (Index i = k - 1; i >= 0; --i) 00172 { 00173 x_new(i) += y(i); 00174 // apply Householder reflection H_{i} to x_new 00175 x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); 00176 } 00177 00178 x += x_new; 00179 00180 if(stop) 00181 { 00182 return true; 00183 } 00184 else 00185 { 00186 k=0; 00187 00188 // reset data for restart 00189 p0.noalias() = rhs - mat*x; 00190 r0 = precond.solve(p0); 00191 00192 // clear Hessenberg matrix and Householder data 00193 H.setZero(); 00194 w.setZero(); 00195 tau.setZero(); 00196 00197 // generate first Householder vector 00198 r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta); 00199 w(0) = Scalar(beta); 00200 } 00201 } 00202 } 00203 00204 return false; 00205 00206 } 00207 00208 } 00209 00210 template< typename _MatrixType, 00211 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> > 00212 class GMRES; 00213 00214 namespace internal { 00215 00216 template< typename _MatrixType, typename _Preconditioner> 00217 struct traits<GMRES<_MatrixType,_Preconditioner> > 00218 { 00219 typedef _MatrixType MatrixType; 00220 typedef _Preconditioner Preconditioner; 00221 }; 00222 00223 } 00224 00259 template< typename _MatrixType, typename _Preconditioner> 00260 class GMRES : public IterativeSolverBase<GMRES<_MatrixType,_Preconditioner> > 00261 { 00262 typedef IterativeSolverBase<GMRES> Base; 00263 using Base::matrix; 00264 using Base::m_error; 00265 using Base::m_iterations; 00266 using Base::m_info; 00267 using Base::m_isInitialized; 00268 00269 private: 00270 Index m_restart; 00271 00272 public: 00273 using Base::_solve_impl; 00274 typedef _MatrixType MatrixType; 00275 typedef typename MatrixType::Scalar Scalar; 00276 typedef typename MatrixType::RealScalar RealScalar; 00277 typedef _Preconditioner Preconditioner; 00278 00279 public: 00280 00282 GMRES() : Base(), m_restart(30) {} 00283 00294 template<typename MatrixDerived> 00295 explicit GMRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_restart(30) {} 00296 00297 ~GMRES() {} 00298 00301 Index get_restart() { return m_restart; } 00302 00306 void set_restart(const Index restart) { m_restart=restart; } 00307 00309 template<typename Rhs,typename Dest> 00310 void _solve_with_guess_impl(const Rhs& b, Dest& x) const 00311 { 00312 bool failed = false; 00313 for(Index j=0; j<b.cols(); ++j) 00314 { 00315 m_iterations = Base::maxIterations(); 00316 m_error = Base::m_tolerance; 00317 00318 typename Dest::ColXpr xj(x,j); 00319 if(!internal::gmres(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error)) 00320 failed = true; 00321 } 00322 m_info = failed ? NumericalIssue 00323 : m_error <= Base::m_tolerance ? Success 00324 : NoConvergence; 00325 m_isInitialized = true; 00326 } 00327 00329 template<typename Rhs,typename Dest> 00330 void _solve_impl(const Rhs& b, MatrixBase<Dest> &x) const 00331 { 00332 x = b; 00333 if(x.squaredNorm() == 0) return; // Check Zero right hand side 00334 _solve_with_guess_impl(b,x.derived()); 00335 } 00336 00337 protected: 00338 00339 }; 00340 00341 } // end namespace Eigen 00342 00343 #endif // EIGEN_GMRES_H