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. 250 * 251 * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner 252 */ 253 template< typename _MatrixType, typename _Preconditioner> 254 class GMRES : public IterativeSolverBase<GMRES<_MatrixType,_Preconditioner> > 255 { 256 typedef IterativeSolverBase<GMRES> Base; 257 using Base::mp_matrix; 258 using Base::m_error; 259 using Base::m_iterations; 260 using Base::m_info; 261 using Base::m_isInitialized; 262 263 private: 264 int m_restart; 265 266 public: 267 typedef _MatrixType MatrixType; 268 typedef typename MatrixType::Scalar Scalar; 269 typedef typename MatrixType::Index Index; 270 typedef typename MatrixType::RealScalar RealScalar; 271 typedef _Preconditioner Preconditioner; 272 273 public: 274 275 /** Default constructor. */ 276 GMRES() : Base(), m_restart(30) {} 277 278 /** Initialize the solver with matrix \a A for further \c Ax=b solving. 279 * 280 * This constructor is a shortcut for the default constructor followed 281 * by a call to compute(). 282 * 283 * \warning this class stores a reference to the matrix A as well as some 284 * precomputed values that depend on it. Therefore, if \a A is changed 285 * this class becomes invalid. Call compute() to update it with the new 286 * matrix A, or modify a copy of A. 287 */ 288 GMRES(const MatrixType& A) : Base(A), m_restart(30) {} 289 290 ~GMRES() {} 291 292 /** Get the number of iterations after that a restart is performed. 293 */ 294 int get_restart() { return m_restart; } 295 296 /** Set the number of iterations after that a restart is performed. 297 * \param restart number of iterations for a restarti, default is 30. 298 */ 299 void set_restart(const int restart) { m_restart=restart; } 300 301 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A 302 * \a x0 as an initial solution. 303 * 304 * \sa compute() 305 */ 306 template<typename Rhs,typename Guess> 307 inline const internal::solve_retval_with_guess<GMRES, Rhs, Guess> 308 solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const 309 { 310 eigen_assert(m_isInitialized && "GMRES is not initialized."); 311 eigen_assert(Base::rows()==b.rows() 312 && "GMRES::solve(): invalid number of rows of the right hand side matrix b"); 313 return internal::solve_retval_with_guess 314 <GMRES, Rhs, Guess>(*this, b.derived(), x0); 315 } 316 317 /** \internal */ 318 template<typename Rhs,typename Dest> 319 void _solveWithGuess(const Rhs& b, Dest& x) const 320 { 321 bool failed = false; 322 for(int j=0; j<b.cols(); ++j) 323 { 324 m_iterations = Base::maxIterations(); 325 m_error = Base::m_tolerance; 326 327 typename Dest::ColXpr xj(x,j); 328 if(!internal::gmres(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error)) 329 failed = true; 330 } 331 m_info = failed ? NumericalIssue 332 : m_error <= Base::m_tolerance ? Success 333 : NoConvergence; 334 m_isInitialized = true; 335 } 336 337 /** \internal */ 338 template<typename Rhs,typename Dest> 339 void _solve(const Rhs& b, Dest& x) const 340 { 341 x = b; 342 if(x.squaredNorm() == 0) return; // Check Zero right hand side 343 _solveWithGuess(b,x); 344 } 345 346 protected: 347 348 }; 349 350 351 namespace internal { 352 353 template<typename _MatrixType, typename _Preconditioner, typename Rhs> 354 struct solve_retval<GMRES<_MatrixType, _Preconditioner>, Rhs> 355 : solve_retval_base<GMRES<_MatrixType, _Preconditioner>, Rhs> 356 { 357 typedef GMRES<_MatrixType, _Preconditioner> Dec; 358 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 359 360 template<typename Dest> void evalTo(Dest& dst) const 361 { 362 dec()._solve(rhs(),dst); 363 } 364 }; 365 366 } // end namespace internal 367 368 } // end namespace Eigen 369 370 #endif // EIGEN_GMRES_H 371