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