1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_DGMRES_H 11 #define EIGEN_DGMRES_H 12 13 #include <Eigen/Eigenvalues> 14 15 namespace Eigen { 16 17 template< typename _MatrixType, 18 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> > 19 class DGMRES; 20 21 namespace internal { 22 23 template< typename _MatrixType, typename _Preconditioner> 24 struct traits<DGMRES<_MatrixType,_Preconditioner> > 25 { 26 typedef _MatrixType MatrixType; 27 typedef _Preconditioner Preconditioner; 28 }; 29 30 /** \brief Computes a permutation vector to have a sorted sequence 31 * \param vec The vector to reorder. 32 * \param perm gives the sorted sequence on output. Must be initialized with 0..n-1 33 * \param ncut Put the ncut smallest elements at the end of the vector 34 * WARNING This is an expensive sort, so should be used only 35 * for small size vectors 36 * TODO Use modified QuickSplit or std::nth_element to get the smallest values 37 */ 38 template <typename VectorType, typename IndexType> 39 void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType::Scalar& ncut) 40 { 41 eigen_assert(vec.size() == perm.size()); 42 typedef typename IndexType::Scalar Index; 43 bool flag; 44 for (Index k = 0; k < ncut; k++) 45 { 46 flag = false; 47 for (Index j = 0; j < vec.size()-1; j++) 48 { 49 if ( vec(perm(j)) < vec(perm(j+1)) ) 50 { 51 std::swap(perm(j),perm(j+1)); 52 flag = true; 53 } 54 if (!flag) break; // The vector is in sorted order 55 } 56 } 57 } 58 59 } 60 /** 61 * \ingroup IterativeLInearSolvers_Module 62 * \brief A Restarted GMRES with deflation. 63 * This class implements a modification of the GMRES solver for 64 * sparse linear systems. The basis is built with modified 65 * Gram-Schmidt. At each restart, a few approximated eigenvectors 66 * corresponding to the smallest eigenvalues are used to build a 67 * preconditioner for the next cycle. This preconditioner 68 * for deflation can be combined with any other preconditioner, 69 * the IncompleteLUT for instance. The preconditioner is applied 70 * at right of the matrix and the combination is multiplicative. 71 * 72 * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix. 73 * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner 74 * Typical usage : 75 * \code 76 * SparseMatrix<double> A; 77 * VectorXd x, b; 78 * //Fill A and b ... 79 * DGMRES<SparseMatrix<double> > solver; 80 * solver.set_restart(30); // Set restarting value 81 * solver.setEigenv(1); // Set the number of eigenvalues to deflate 82 * solver.compute(A); 83 * x = solver.solve(b); 84 * \endcode 85 * 86 * DGMRES can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink. 87 * 88 * References : 89 * [1] D. NUENTSA WAKAM and F. PACULL, Memory Efficient Hybrid 90 * Algebraic Solvers for Linear Systems Arising from Compressible 91 * Flows, Computers and Fluids, In Press, 92 * http://dx.doi.org/10.1016/j.compfluid.2012.03.023 93 * [2] K. Burrage and J. Erhel, On the performance of various 94 * adaptive preconditioned GMRES strategies, 5(1998), 101-121. 95 * [3] J. Erhel, K. Burrage and B. Pohl, Restarted GMRES 96 * preconditioned by deflation,J. Computational and Applied 97 * Mathematics, 69(1996), 303-318. 98 99 * 100 */ 101 template< typename _MatrixType, typename _Preconditioner> 102 class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > 103 { 104 typedef IterativeSolverBase<DGMRES> Base; 105 using Base::matrix; 106 using Base::m_error; 107 using Base::m_iterations; 108 using Base::m_info; 109 using Base::m_isInitialized; 110 using Base::m_tolerance; 111 public: 112 using Base::_solve_impl; 113 typedef _MatrixType MatrixType; 114 typedef typename MatrixType::Scalar Scalar; 115 typedef typename MatrixType::Index Index; 116 typedef typename MatrixType::StorageIndex StorageIndex; 117 typedef typename MatrixType::RealScalar RealScalar; 118 typedef _Preconditioner Preconditioner; 119 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 120 typedef Matrix<RealScalar,Dynamic,Dynamic> DenseRealMatrix; 121 typedef Matrix<Scalar,Dynamic,1> DenseVector; 122 typedef Matrix<RealScalar,Dynamic,1> DenseRealVector; 123 typedef Matrix<std::complex<RealScalar>, Dynamic, 1> ComplexVector; 124 125 126 /** Default constructor. */ 127 DGMRES() : Base(),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {} 128 129 /** Initialize the solver with matrix \a A for further \c Ax=b solving. 130 * 131 * This constructor is a shortcut for the default constructor followed 132 * by a call to compute(). 133 * 134 * \warning this class stores a reference to the matrix A as well as some 135 * precomputed values that depend on it. Therefore, if \a A is changed 136 * this class becomes invalid. Call compute() to update it with the new 137 * matrix A, or modify a copy of A. 138 */ 139 template<typename MatrixDerived> 140 explicit DGMRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {} 141 142 ~DGMRES() {} 143 144 /** \internal */ 145 template<typename Rhs,typename Dest> 146 void _solve_with_guess_impl(const Rhs& b, Dest& x) const 147 { 148 bool failed = false; 149 for(int j=0; j<b.cols(); ++j) 150 { 151 m_iterations = Base::maxIterations(); 152 m_error = Base::m_tolerance; 153 154 typename Dest::ColXpr xj(x,j); 155 dgmres(matrix(), b.col(j), xj, Base::m_preconditioner); 156 } 157 m_info = failed ? NumericalIssue 158 : m_error <= Base::m_tolerance ? Success 159 : NoConvergence; 160 m_isInitialized = true; 161 } 162 163 /** \internal */ 164 template<typename Rhs,typename Dest> 165 void _solve_impl(const Rhs& b, MatrixBase<Dest>& x) const 166 { 167 x = b; 168 _solve_with_guess_impl(b,x.derived()); 169 } 170 /** 171 * Get the restart value 172 */ 173 int restart() { return m_restart; } 174 175 /** 176 * Set the restart value (default is 30) 177 */ 178 void set_restart(const int restart) { m_restart=restart; } 179 180 /** 181 * Set the number of eigenvalues to deflate at each restart 182 */ 183 void setEigenv(const int neig) 184 { 185 m_neig = neig; 186 if (neig+1 > m_maxNeig) m_maxNeig = neig+1; // To allow for complex conjugates 187 } 188 189 /** 190 * Get the size of the deflation subspace size 191 */ 192 int deflSize() {return m_r; } 193 194 /** 195 * Set the maximum size of the deflation subspace 196 */ 197 void setMaxEigenv(const int maxNeig) { m_maxNeig = maxNeig; } 198 199 protected: 200 // DGMRES algorithm 201 template<typename Rhs, typename Dest> 202 void dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, const Preconditioner& precond) const; 203 // Perform one cycle of GMRES 204 template<typename Dest> 205 int dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const; 206 // Compute data to use for deflation 207 int dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const; 208 // Apply deflation to a vector 209 template<typename RhsType, typename DestType> 210 int dgmresApplyDeflation(const RhsType& In, DestType& Out) const; 211 ComplexVector schurValues(const ComplexSchur<DenseMatrix>& schurofH) const; 212 ComplexVector schurValues(const RealSchur<DenseMatrix>& schurofH) const; 213 // Init data for deflation 214 void dgmresInitDeflation(Index& rows) const; 215 mutable DenseMatrix m_V; // Krylov basis vectors 216 mutable DenseMatrix m_H; // Hessenberg matrix 217 mutable DenseMatrix m_Hes; // Initial hessenberg matrix wihout Givens rotations applied 218 mutable Index m_restart; // Maximum size of the Krylov subspace 219 mutable DenseMatrix m_U; // Vectors that form the basis of the invariant subspace 220 mutable DenseMatrix m_MU; // matrix operator applied to m_U (for next cycles) 221 mutable DenseMatrix m_T; /* T=U^T*M^{-1}*A*U */ 222 mutable PartialPivLU<DenseMatrix> m_luT; // LU factorization of m_T 223 mutable StorageIndex m_neig; //Number of eigenvalues to extract at each restart 224 mutable int m_r; // Current number of deflated eigenvalues, size of m_U 225 mutable int m_maxNeig; // Maximum number of eigenvalues to deflate 226 mutable RealScalar m_lambdaN; //Modulus of the largest eigenvalue of A 227 mutable bool m_isDeflAllocated; 228 mutable bool m_isDeflInitialized; 229 230 //Adaptive strategy 231 mutable RealScalar m_smv; // Smaller multiple of the remaining number of steps allowed 232 mutable bool m_force; // Force the use of deflation at each restart 233 234 }; 235 /** 236 * \brief Perform several cycles of restarted GMRES with modified Gram Schmidt, 237 * 238 * A right preconditioner is used combined with deflation. 239 * 240 */ 241 template< typename _MatrixType, typename _Preconditioner> 242 template<typename Rhs, typename Dest> 243 void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, 244 const Preconditioner& precond) const 245 { 246 //Initialization 247 int n = mat.rows(); 248 DenseVector r0(n); 249 int nbIts = 0; 250 m_H.resize(m_restart+1, m_restart); 251 m_Hes.resize(m_restart, m_restart); 252 m_V.resize(n,m_restart+1); 253 //Initial residual vector and intial norm 254 x = precond.solve(x); 255 r0 = rhs - mat * x; 256 RealScalar beta = r0.norm(); 257 RealScalar normRhs = rhs.norm(); 258 m_error = beta/normRhs; 259 if(m_error < m_tolerance) 260 m_info = Success; 261 else 262 m_info = NoConvergence; 263 264 // Iterative process 265 while (nbIts < m_iterations && m_info == NoConvergence) 266 { 267 dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts); 268 269 // Compute the new residual vector for the restart 270 if (nbIts < m_iterations && m_info == NoConvergence) 271 r0 = rhs - mat * x; 272 } 273 } 274 275 /** 276 * \brief Perform one restart cycle of DGMRES 277 * \param mat The coefficient matrix 278 * \param precond The preconditioner 279 * \param x the new approximated solution 280 * \param r0 The initial residual vector 281 * \param beta The norm of the residual computed so far 282 * \param normRhs The norm of the right hand side vector 283 * \param nbIts The number of iterations 284 */ 285 template< typename _MatrixType, typename _Preconditioner> 286 template<typename Dest> 287 int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const 288 { 289 //Initialization 290 DenseVector g(m_restart+1); // Right hand side of the least square problem 291 g.setZero(); 292 g(0) = Scalar(beta); 293 m_V.col(0) = r0/beta; 294 m_info = NoConvergence; 295 std::vector<JacobiRotation<Scalar> >gr(m_restart); // Givens rotations 296 int it = 0; // Number of inner iterations 297 int n = mat.rows(); 298 DenseVector tv1(n), tv2(n); //Temporary vectors 299 while (m_info == NoConvergence && it < m_restart && nbIts < m_iterations) 300 { 301 // Apply preconditioner(s) at right 302 if (m_isDeflInitialized ) 303 { 304 dgmresApplyDeflation(m_V.col(it), tv1); // Deflation 305 tv2 = precond.solve(tv1); 306 } 307 else 308 { 309 tv2 = precond.solve(m_V.col(it)); // User's selected preconditioner 310 } 311 tv1 = mat * tv2; 312 313 // Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt 314 Scalar coef; 315 for (int i = 0; i <= it; ++i) 316 { 317 coef = tv1.dot(m_V.col(i)); 318 tv1 = tv1 - coef * m_V.col(i); 319 m_H(i,it) = coef; 320 m_Hes(i,it) = coef; 321 } 322 // Normalize the vector 323 coef = tv1.norm(); 324 m_V.col(it+1) = tv1/coef; 325 m_H(it+1, it) = coef; 326 // m_Hes(it+1,it) = coef; 327 328 // FIXME Check for happy breakdown 329 330 // Update Hessenberg matrix with Givens rotations 331 for (int i = 1; i <= it; ++i) 332 { 333 m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint()); 334 } 335 // Compute the new plane rotation 336 gr[it].makeGivens(m_H(it, it), m_H(it+1,it)); 337 // Apply the new rotation 338 m_H.col(it).applyOnTheLeft(it,it+1,gr[it].adjoint()); 339 g.applyOnTheLeft(it,it+1, gr[it].adjoint()); 340 341 beta = std::abs(g(it+1)); 342 m_error = beta/normRhs; 343 // std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl; 344 it++; nbIts++; 345 346 if (m_error < m_tolerance) 347 { 348 // The method has converged 349 m_info = Success; 350 break; 351 } 352 } 353 354 // Compute the new coefficients by solving the least square problem 355 // it++; 356 //FIXME Check first if the matrix is singular ... zero diagonal 357 DenseVector nrs(m_restart); 358 nrs = m_H.topLeftCorner(it,it).template triangularView<Upper>().solve(g.head(it)); 359 360 // Form the new solution 361 if (m_isDeflInitialized) 362 { 363 tv1 = m_V.leftCols(it) * nrs; 364 dgmresApplyDeflation(tv1, tv2); 365 x = x + precond.solve(tv2); 366 } 367 else 368 x = x + precond.solve(m_V.leftCols(it) * nrs); 369 370 // Go for a new cycle and compute data for deflation 371 if(nbIts < m_iterations && m_info == NoConvergence && m_neig > 0 && (m_r+m_neig) < m_maxNeig) 372 dgmresComputeDeflationData(mat, precond, it, m_neig); 373 return 0; 374 375 } 376 377 378 template< typename _MatrixType, typename _Preconditioner> 379 void DGMRES<_MatrixType, _Preconditioner>::dgmresInitDeflation(Index& rows) const 380 { 381 m_U.resize(rows, m_maxNeig); 382 m_MU.resize(rows, m_maxNeig); 383 m_T.resize(m_maxNeig, m_maxNeig); 384 m_lambdaN = 0.0; 385 m_isDeflAllocated = true; 386 } 387 388 template< typename _MatrixType, typename _Preconditioner> 389 inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const ComplexSchur<DenseMatrix>& schurofH) const 390 { 391 return schurofH.matrixT().diagonal(); 392 } 393 394 template< typename _MatrixType, typename _Preconditioner> 395 inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const RealSchur<DenseMatrix>& schurofH) const 396 { 397 typedef typename MatrixType::Index Index; 398 const DenseMatrix& T = schurofH.matrixT(); 399 Index it = T.rows(); 400 ComplexVector eig(it); 401 Index j = 0; 402 while (j < it-1) 403 { 404 if (T(j+1,j) ==Scalar(0)) 405 { 406 eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0)); 407 j++; 408 } 409 else 410 { 411 eig(j) = std::complex<RealScalar>(T(j,j),T(j+1,j)); 412 eig(j+1) = std::complex<RealScalar>(T(j,j+1),T(j+1,j+1)); 413 j++; 414 } 415 } 416 if (j < it-1) eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0)); 417 return eig; 418 } 419 420 template< typename _MatrixType, typename _Preconditioner> 421 int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const 422 { 423 // First, find the Schur form of the Hessenberg matrix H 424 typename internal::conditional<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> >::type schurofH; 425 bool computeU = true; 426 DenseMatrix matrixQ(it,it); 427 matrixQ.setIdentity(); 428 schurofH.computeFromHessenberg(m_Hes.topLeftCorner(it,it), matrixQ, computeU); 429 430 ComplexVector eig(it); 431 Matrix<StorageIndex,Dynamic,1>perm(it); 432 eig = this->schurValues(schurofH); 433 434 // Reorder the absolute values of Schur values 435 DenseRealVector modulEig(it); 436 for (int j=0; j<it; ++j) modulEig(j) = std::abs(eig(j)); 437 perm.setLinSpaced(it,0,it-1); 438 internal::sortWithPermutation(modulEig, perm, neig); 439 440 if (!m_lambdaN) 441 { 442 m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN); 443 } 444 //Count the real number of extracted eigenvalues (with complex conjugates) 445 int nbrEig = 0; 446 while (nbrEig < neig) 447 { 448 if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++; 449 else nbrEig += 2; 450 } 451 // Extract the Schur vectors corresponding to the smallest Ritz values 452 DenseMatrix Sr(it, nbrEig); 453 Sr.setZero(); 454 for (int j = 0; j < nbrEig; j++) 455 { 456 Sr.col(j) = schurofH.matrixU().col(perm(it-j-1)); 457 } 458 459 // Form the Schur vectors of the initial matrix using the Krylov basis 460 DenseMatrix X; 461 X = m_V.leftCols(it) * Sr; 462 if (m_r) 463 { 464 // Orthogonalize X against m_U using modified Gram-Schmidt 465 for (int j = 0; j < nbrEig; j++) 466 for (int k =0; k < m_r; k++) 467 X.col(j) = X.col(j) - (m_U.col(k).dot(X.col(j)))*m_U.col(k); 468 } 469 470 // Compute m_MX = A * M^-1 * X 471 Index m = m_V.rows(); 472 if (!m_isDeflAllocated) 473 dgmresInitDeflation(m); 474 DenseMatrix MX(m, nbrEig); 475 DenseVector tv1(m); 476 for (int j = 0; j < nbrEig; j++) 477 { 478 tv1 = mat * X.col(j); 479 MX.col(j) = precond.solve(tv1); 480 } 481 482 //Update m_T = [U'MU U'MX; X'MU X'MX] 483 m_T.block(m_r, m_r, nbrEig, nbrEig) = X.transpose() * MX; 484 if(m_r) 485 { 486 m_T.block(0, m_r, m_r, nbrEig) = m_U.leftCols(m_r).transpose() * MX; 487 m_T.block(m_r, 0, nbrEig, m_r) = X.transpose() * m_MU.leftCols(m_r); 488 } 489 490 // Save X into m_U and m_MX in m_MU 491 for (int j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j); 492 for (int j = 0; j < nbrEig; j++) m_MU.col(m_r+j) = MX.col(j); 493 // Increase the size of the invariant subspace 494 m_r += nbrEig; 495 496 // Factorize m_T into m_luT 497 m_luT.compute(m_T.topLeftCorner(m_r, m_r)); 498 499 //FIXME CHeck if the factorization was correctly done (nonsingular matrix) 500 m_isDeflInitialized = true; 501 return 0; 502 } 503 template<typename _MatrixType, typename _Preconditioner> 504 template<typename RhsType, typename DestType> 505 int DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x, DestType &y) const 506 { 507 DenseVector x1 = m_U.leftCols(m_r).transpose() * x; 508 y = x + m_U.leftCols(m_r) * ( m_lambdaN * m_luT.solve(x1) - x1); 509 return 0; 510 } 511 512 } // end namespace Eigen 513 #endif 514