1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // Copyright (C) 2012, 2014 Kolja Brix <brix (at) igpm.rwth-aaachen.de> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_GMRES_H 12 #define EIGEN_GMRES_H 13 14 namespace Eigen { 15 16 namespace internal { 17 18 /** 19 * Generalized Minimal Residual Algorithm based on the 20 * Arnoldi algorithm implemented with Householder reflections. 21 * 22 * Parameters: 23 * \param mat matrix of linear system of equations 24 * \param Rhs right hand side vector of linear system of equations 25 * \param x on input: initial guess, on output: solution 26 * \param precond preconditioner used 27 * \param iters on input: maximum number of iterations to perform 28 * on output: number of iterations performed 29 * \param restart number of iterations for a restart 30 * \param tol_error on input: residual tolerance 31 * on output: residuum achieved 32 * 33 * \sa IterativeMethods::bicgstab() 34 * 35 * 36 * For references, please see: 37 * 38 * Saad, Y. and Schultz, M. H. 39 * GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. 40 * SIAM J.Sci.Stat.Comp. 7, 1986, pp. 856 - 869. 41 * 42 * Saad, Y. 43 * Iterative Methods for Sparse Linear Systems. 44 * Society for Industrial and Applied Mathematics, Philadelphia, 2003. 45 * 46 * Walker, H. F. 47 * Implementations of the GMRES method. 48 * Comput.Phys.Comm. 53, 1989, pp. 311 - 320. 49 * 50 * Walker, H. F. 51 * Implementation of the GMRES Method using Householder Transformations. 52 * SIAM J.Sci.Stat.Comp. 9, 1988, pp. 152 - 163. 53 * 54 */ 55 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner> 56 bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Preconditioner & precond, 57 int &iters, const int &restart, typename Dest::RealScalar & tol_error) { 58 59 using std::sqrt; 60 using std::abs; 61 62 typedef typename Dest::RealScalar RealScalar; 63 typedef typename Dest::Scalar Scalar; 64 typedef Matrix < Scalar, Dynamic, 1 > VectorType; 65 typedef Matrix < Scalar, Dynamic, Dynamic > FMatrixType; 66 67 RealScalar tol = tol_error; 68 const int maxIters = iters; 69 iters = 0; 70 71 const int m = mat.rows(); 72 73 VectorType p0 = rhs - mat*x; 74 VectorType r0 = precond.solve(p0); 75 76 // is initial guess already good enough? 77 if(abs(r0.norm()) < tol) { 78 return true; 79 } 80 81 VectorType w = VectorType::Zero(restart + 1); 82 83 FMatrixType H = FMatrixType::Zero(m, restart + 1); // Hessenberg matrix 84 VectorType tau = VectorType::Zero(restart + 1); 85 std::vector < JacobiRotation < Scalar > > G(restart); 86 87 // generate first Householder vector 88 VectorType e(m-1); 89 RealScalar beta; 90 r0.makeHouseholder(e, tau.coeffRef(0), beta); 91 w(0)=(Scalar) beta; 92 H.bottomLeftCorner(m - 1, 1) = e; 93 94 for (int k = 1; k <= restart; ++k) { 95 96 ++iters; 97 98 VectorType v = VectorType::Unit(m, k - 1), workspace(m); 99 100 // apply Householder reflections H_{1} ... H_{k-1} to v 101 for (int i = k - 1; i >= 0; --i) { 102 v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); 103 } 104 105 // apply matrix M to v: v = mat * v; 106 VectorType t=mat*v; 107 v=precond.solve(t); 108 109 // apply Householder reflections H_{k-1} ... H_{1} to v 110 for (int i = 0; i < k; ++i) { 111 v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); 112 } 113 114 if (v.tail(m - k).norm() != 0.0) { 115 116 if (k <= restart) { 117 118 // generate new Householder vector 119 VectorType e(m - k - 1); 120 RealScalar beta; 121 v.tail(m - k).makeHouseholder(e, tau.coeffRef(k), beta); 122 H.col(k).tail(m - k - 1) = e; 123 124 // apply Householder reflection H_{k} to v 125 v.tail(m - k).applyHouseholderOnTheLeft(H.col(k).tail(m - k - 1), tau.coeffRef(k), workspace.data()); 126 127 } 128 } 129 130 if (k > 1) { 131 for (int i = 0; i < k - 1; ++i) { 132 // apply old Givens rotations to v 133 v.applyOnTheLeft(i, i + 1, G[i].adjoint()); 134 } 135 } 136 137 if (k<m && v(k) != (Scalar) 0) { 138 // determine next Givens rotation 139 G[k - 1].makeGivens(v(k - 1), v(k)); 140 141 // apply Givens rotation to v and w 142 v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint()); 143 w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint()); 144 145 } 146 147 // insert coefficients into upper matrix triangle 148 H.col(k - 1).head(k) = v.head(k); 149 150 bool stop=(k==m || abs(w(k)) < tol || iters == maxIters); 151 152 if (stop || k == restart) { 153 154 // solve upper triangular system 155 VectorType y = w.head(k); 156 H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y); 157 158 // use Horner-like scheme to calculate solution vector 159 VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1); 160 161 // apply Householder reflection H_{k} to x_new 162 x_new.tail(m - k + 1).applyHouseholderOnTheLeft(H.col(k - 1).tail(m - k), tau.coeffRef(k - 1), workspace.data()); 163 164 for (int i = k - 2; i >= 0; --i) { 165 x_new += y(i) * VectorType::Unit(m, i); 166 // apply Householder reflection H_{i} to x_new 167 x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); 168 } 169 170 x += x_new; 171 172 if (stop) { 173 return true; 174 } else { 175 k=0; 176 177 // reset data for a restart r0 = rhs - mat * x; 178 VectorType p0=mat*x; 179 VectorType p1=precond.solve(p0); 180 r0 = rhs - p1; 181 // r0_sqnorm = r0.squaredNorm(); 182 w = VectorType::Zero(restart + 1); 183 H = FMatrixType::Zero(m, restart + 1); 184 tau = VectorType::Zero(restart + 1); 185 186 // generate first Householder vector 187 RealScalar beta; 188 r0.makeHouseholder(e, tau.coeffRef(0), beta); 189 w(0)=(Scalar) beta; 190 H.bottomLeftCorner(m - 1, 1) = e; 191 192 } 193 194 } 195 196 197 198 } 199 200 return false; 201 202 } 203 204 } 205 206 template< typename _MatrixType, 207 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> > 208 class GMRES; 209 210 namespace internal { 211 212 template< typename _MatrixType, typename _Preconditioner> 213 struct traits<GMRES<_MatrixType,_Preconditioner> > 214 { 215 typedef _MatrixType MatrixType; 216 typedef _Preconditioner Preconditioner; 217 }; 218 219 } 220 221 /** \ingroup IterativeLinearSolvers_Module 222 * \brief A GMRES solver for sparse square problems 223 * 224 * This class allows to solve for A.x = b sparse linear problems using a generalized minimal 225 * residual method. The vectors x and b can be either dense or sparse. 226 * 227 * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix. 228 * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner 229 * 230 * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations() 231 * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations 232 * and NumTraits<Scalar>::epsilon() for the tolerance. 233 * 234 * This class can be used as the direct solver classes. Here is a typical usage example: 235 * \code 236 * int n = 10000; 237 * VectorXd x(n), b(n); 238 * SparseMatrix<double> A(n,n); 239 * // fill A and b 240 * GMRES<SparseMatrix<double> > solver(A); 241 * x = solver.solve(b); 242 * std::cout << "#iterations: " << solver.iterations() << std::endl; 243 * std::cout << "estimated error: " << solver.error() << std::endl; 244 * // update b, and solve again 245 * x = solver.solve(b); 246 * \endcode 247 * 248 * By default the iterations start with x=0 as an initial guess of the solution. 249 * One can control the start using the solveWithGuess() method. Here is a step by 250 * step execution example starting with a random guess and printing the evolution 251 * of the estimated error: 252 * * \code 253 * x = VectorXd::Random(n); 254 * solver.setMaxIterations(1); 255 * int i = 0; 256 * do { 257 * x = solver.solveWithGuess(b,x); 258 * std::cout << i << " : " << solver.error() << std::endl; 259 * ++i; 260 * } while (solver.info()!=Success && i<100); 261 * \endcode 262 * Note that such a step by step excution is slightly slower. 263 * 264 * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner 265 */ 266 template< typename _MatrixType, typename _Preconditioner> 267 class GMRES : public IterativeSolverBase<GMRES<_MatrixType,_Preconditioner> > 268 { 269 typedef IterativeSolverBase<GMRES> Base; 270 using Base::mp_matrix; 271 using Base::m_error; 272 using Base::m_iterations; 273 using Base::m_info; 274 using Base::m_isInitialized; 275 276 private: 277 int m_restart; 278 279 public: 280 typedef _MatrixType MatrixType; 281 typedef typename MatrixType::Scalar Scalar; 282 typedef typename MatrixType::Index Index; 283 typedef typename MatrixType::RealScalar RealScalar; 284 typedef _Preconditioner Preconditioner; 285 286 public: 287 288 /** Default constructor. */ 289 GMRES() : Base(), m_restart(30) {} 290 291 /** Initialize the solver with matrix \a A for further \c Ax=b solving. 292 * 293 * This constructor is a shortcut for the default constructor followed 294 * by a call to compute(). 295 * 296 * \warning this class stores a reference to the matrix A as well as some 297 * precomputed values that depend on it. Therefore, if \a A is changed 298 * this class becomes invalid. Call compute() to update it with the new 299 * matrix A, or modify a copy of A. 300 */ 301 GMRES(const MatrixType& A) : Base(A), m_restart(30) {} 302 303 ~GMRES() {} 304 305 /** Get the number of iterations after that a restart is performed. 306 */ 307 int get_restart() { return m_restart; } 308 309 /** Set the number of iterations after that a restart is performed. 310 * \param restart number of iterations for a restarti, default is 30. 311 */ 312 void set_restart(const int restart) { m_restart=restart; } 313 314 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A 315 * \a x0 as an initial solution. 316 * 317 * \sa compute() 318 */ 319 template<typename Rhs,typename Guess> 320 inline const internal::solve_retval_with_guess<GMRES, Rhs, Guess> 321 solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const 322 { 323 eigen_assert(m_isInitialized && "GMRES is not initialized."); 324 eigen_assert(Base::rows()==b.rows() 325 && "GMRES::solve(): invalid number of rows of the right hand side matrix b"); 326 return internal::solve_retval_with_guess 327 <GMRES, Rhs, Guess>(*this, b.derived(), x0); 328 } 329 330 /** \internal */ 331 template<typename Rhs,typename Dest> 332 void _solveWithGuess(const Rhs& b, Dest& x) const 333 { 334 bool failed = false; 335 for(int j=0; j<b.cols(); ++j) 336 { 337 m_iterations = Base::maxIterations(); 338 m_error = Base::m_tolerance; 339 340 typename Dest::ColXpr xj(x,j); 341 if(!internal::gmres(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error)) 342 failed = true; 343 } 344 m_info = failed ? NumericalIssue 345 : m_error <= Base::m_tolerance ? Success 346 : NoConvergence; 347 m_isInitialized = true; 348 } 349 350 /** \internal */ 351 template<typename Rhs,typename Dest> 352 void _solve(const Rhs& b, Dest& x) const 353 { 354 x = b; 355 if(x.squaredNorm() == 0) return; // Check Zero right hand side 356 _solveWithGuess(b,x); 357 } 358 359 protected: 360 361 }; 362 363 364 namespace internal { 365 366 template<typename _MatrixType, typename _Preconditioner, typename Rhs> 367 struct solve_retval<GMRES<_MatrixType, _Preconditioner>, Rhs> 368 : solve_retval_base<GMRES<_MatrixType, _Preconditioner>, Rhs> 369 { 370 typedef GMRES<_MatrixType, _Preconditioner> Dec; 371 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 372 373 template<typename Dest> void evalTo(Dest& dst) const 374 { 375 dec()._solve(rhs(),dst); 376 } 377 }; 378 379 } // end namespace internal 380 381 } // end namespace Eigen 382 383 #endif // EIGEN_GMRES_H 384
