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