Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /include/eigen3/unsupported/Eigen/src/IterativeSolvers/DGMRES.h was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009 
0010 #ifndef EIGEN_DGMRES_H
0011 #define EIGEN_DGMRES_H
0012 
0013 #include "../../../../Eigen/Eigenvalues"
0014 
0015 namespace Eigen { 
0016   
0017 template< typename _MatrixType,
0018           typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
0019 class DGMRES;
0020 
0021 namespace internal {
0022 
0023 template< typename _MatrixType, typename _Preconditioner>
0024 struct traits<DGMRES<_MatrixType,_Preconditioner> >
0025 {
0026   typedef _MatrixType MatrixType;
0027   typedef _Preconditioner Preconditioner;
0028 };
0029 
0030 /** \brief Computes a permutation vector to have a sorted sequence
0031   * \param vec The vector to reorder.
0032   * \param perm gives the sorted sequence on output. Must be initialized with 0..n-1
0033   * \param ncut Put  the ncut smallest elements at the end of the vector
0034   * WARNING This is an expensive sort, so should be used only 
0035   * for small size vectors
0036   * TODO Use modified QuickSplit or std::nth_element to get the smallest values 
0037   */
0038 template <typename VectorType, typename IndexType>
0039 void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType::Scalar& ncut)
0040 {
0041   eigen_assert(vec.size() == perm.size());
0042   bool flag; 
0043   for (Index k  = 0; k < ncut; k++)
0044   {
0045     flag = false;
0046     for (Index j = 0; j < vec.size()-1; j++)
0047     {
0048       if ( vec(perm(j)) < vec(perm(j+1)) )
0049       {
0050         std::swap(perm(j),perm(j+1)); 
0051         flag = true;
0052       }
0053       if (!flag) break; // The vector is in sorted order
0054     }
0055   }
0056 }
0057 
0058 }
0059 /**
0060  * \ingroup IterativeLinearSolvers_Module
0061  * \brief A Restarted GMRES with deflation.
0062  * This class implements a modification of the GMRES solver for
0063  * sparse linear systems. The basis is built with modified 
0064  * Gram-Schmidt. At each restart, a few approximated eigenvectors
0065  * corresponding to the smallest eigenvalues are used to build a
0066  * preconditioner for the next cycle. This preconditioner 
0067  * for deflation can be combined with any other preconditioner, 
0068  * the IncompleteLUT for instance. The preconditioner is applied 
0069  * at right of the matrix and the combination is multiplicative.
0070  * 
0071  * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix.
0072  * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner
0073  * Typical usage :
0074  * \code
0075  * SparseMatrix<double> A;
0076  * VectorXd x, b; 
0077  * //Fill A and b ...
0078  * DGMRES<SparseMatrix<double> > solver;
0079  * solver.set_restart(30); // Set restarting value
0080  * solver.setEigenv(1); // Set the number of eigenvalues to deflate
0081  * solver.compute(A);
0082  * x = solver.solve(b);
0083  * \endcode
0084  * 
0085  * DGMRES can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink.
0086  *
0087  * References :
0088  * [1] D. NUENTSA WAKAM and F. PACULL, Memory Efficient Hybrid
0089  *  Algebraic Solvers for Linear Systems Arising from Compressible
0090  *  Flows, Computers and Fluids, In Press,
0091  *  https://doi.org/10.1016/j.compfluid.2012.03.023   
0092  * [2] K. Burrage and J. Erhel, On the performance of various 
0093  * adaptive preconditioned GMRES strategies, 5(1998), 101-121.
0094  * [3] J. Erhel, K. Burrage and B. Pohl, Restarted GMRES 
0095  *  preconditioned by deflation,J. Computational and Applied
0096  *  Mathematics, 69(1996), 303-318. 
0097 
0098  * 
0099  */
0100 template< typename _MatrixType, typename _Preconditioner>
0101 class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
0102 {
0103     typedef IterativeSolverBase<DGMRES> Base;
0104     using Base::matrix;
0105     using Base::m_error;
0106     using Base::m_iterations;
0107     using Base::m_info;
0108     using Base::m_isInitialized;
0109     using Base::m_tolerance; 
0110   public:
0111     using Base::_solve_impl;
0112     using Base::_solve_with_guess_impl;
0113     typedef _MatrixType MatrixType;
0114     typedef typename MatrixType::Scalar Scalar;
0115     typedef typename MatrixType::StorageIndex StorageIndex;
0116     typedef typename MatrixType::RealScalar RealScalar;
0117     typedef _Preconditioner Preconditioner;
0118     typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 
0119     typedef Matrix<RealScalar,Dynamic,Dynamic> DenseRealMatrix; 
0120     typedef Matrix<Scalar,Dynamic,1> DenseVector;
0121     typedef Matrix<RealScalar,Dynamic,1> DenseRealVector; 
0122     typedef Matrix<std::complex<RealScalar>, Dynamic, 1> ComplexVector;
0123  
0124     
0125   /** Default constructor. */
0126   DGMRES() : Base(),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {}
0127 
0128   /** Initialize the solver with matrix \a A for further \c Ax=b solving.
0129     * 
0130     * This constructor is a shortcut for the default constructor followed
0131     * by a call to compute().
0132     * 
0133     * \warning this class stores a reference to the matrix A as well as some
0134     * precomputed values that depend on it. Therefore, if \a A is changed
0135     * this class becomes invalid. Call compute() to update it with the new
0136     * matrix A, or modify a copy of A.
0137     */
0138   template<typename MatrixDerived>
0139   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) {}
0140 
0141   ~DGMRES() {}
0142   
0143   /** \internal */
0144   template<typename Rhs,typename Dest>
0145   void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
0146   {
0147     EIGEN_STATIC_ASSERT(Rhs::ColsAtCompileTime==1 || Dest::ColsAtCompileTime==1, YOU_TRIED_CALLING_A_VECTOR_METHOD_ON_A_MATRIX);
0148     
0149     m_iterations = Base::maxIterations();
0150     m_error = Base::m_tolerance;
0151     
0152     dgmres(matrix(), b, x, Base::m_preconditioner);
0153   }
0154 
0155   /** 
0156    * Get the restart value
0157     */
0158   Index restart() { return m_restart; }
0159   
0160   /** 
0161    * Set the restart value (default is 30)  
0162    */
0163   void set_restart(const Index restart) { m_restart=restart; }
0164   
0165   /** 
0166    * Set the number of eigenvalues to deflate at each restart 
0167    */
0168   void setEigenv(const Index neig) 
0169   {
0170     m_neig = neig;
0171     if (neig+1 > m_maxNeig) m_maxNeig = neig+1; // To allow for complex conjugates
0172   }
0173   
0174   /** 
0175    * Get the size of the deflation subspace size
0176    */ 
0177   Index deflSize() {return m_r; }
0178   
0179   /**
0180    * Set the maximum size of the deflation subspace
0181    */
0182   void setMaxEigenv(const Index maxNeig) { m_maxNeig = maxNeig; }
0183   
0184   protected:
0185     // DGMRES algorithm 
0186     template<typename Rhs, typename Dest>
0187     void dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, const Preconditioner& precond) const;
0188     // Perform one cycle of GMRES
0189     template<typename Dest>
0190     Index dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, Index& nbIts) const; 
0191     // Compute data to use for deflation 
0192     Index dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const;
0193     // Apply deflation to a vector
0194     template<typename RhsType, typename DestType>
0195     Index dgmresApplyDeflation(const RhsType& In, DestType& Out) const; 
0196     ComplexVector schurValues(const ComplexSchur<DenseMatrix>& schurofH) const;
0197     ComplexVector schurValues(const RealSchur<DenseMatrix>& schurofH) const;
0198     // Init data for deflation
0199     void dgmresInitDeflation(Index& rows) const; 
0200     mutable DenseMatrix m_V; // Krylov basis vectors
0201     mutable DenseMatrix m_H; // Hessenberg matrix 
0202     mutable DenseMatrix m_Hes; // Initial hessenberg matrix without Givens rotations applied
0203     mutable Index m_restart; // Maximum size of the Krylov subspace
0204     mutable DenseMatrix m_U; // Vectors that form the basis of the invariant subspace 
0205     mutable DenseMatrix m_MU; // matrix operator applied to m_U (for next cycles)
0206     mutable DenseMatrix m_T; /* T=U^T*M^{-1}*A*U */
0207     mutable PartialPivLU<DenseMatrix> m_luT; // LU factorization of m_T
0208     mutable StorageIndex m_neig; //Number of eigenvalues to extract at each restart
0209     mutable Index m_r; // Current number of deflated eigenvalues, size of m_U
0210     mutable Index m_maxNeig; // Maximum number of eigenvalues to deflate
0211     mutable RealScalar m_lambdaN; //Modulus of the largest eigenvalue of A
0212     mutable bool m_isDeflAllocated;
0213     mutable bool m_isDeflInitialized;
0214     
0215     //Adaptive strategy 
0216     mutable RealScalar m_smv; // Smaller multiple of the remaining number of steps allowed
0217     mutable bool m_force; // Force the use of deflation at each restart
0218     
0219 }; 
0220 /** 
0221  * \brief Perform several cycles of restarted GMRES with modified Gram Schmidt, 
0222  * 
0223  * A right preconditioner is used combined with deflation.
0224  * 
0225  */
0226 template< typename _MatrixType, typename _Preconditioner>
0227 template<typename Rhs, typename Dest>
0228 void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x,
0229               const Preconditioner& precond) const
0230 {
0231   const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
0232 
0233   RealScalar normRhs = rhs.norm();
0234   if(normRhs <= considerAsZero) 
0235   {
0236     x.setZero();
0237     m_error = 0;
0238     return;
0239   }
0240 
0241   //Initialization
0242   m_isDeflInitialized = false;
0243   Index n = mat.rows(); 
0244   DenseVector r0(n); 
0245   Index nbIts = 0; 
0246   m_H.resize(m_restart+1, m_restart);
0247   m_Hes.resize(m_restart, m_restart);
0248   m_V.resize(n,m_restart+1);
0249   //Initial residual vector and initial norm
0250   if(x.squaredNorm()==0) 
0251     x = precond.solve(rhs);
0252   r0 = rhs - mat * x; 
0253   RealScalar beta = r0.norm(); 
0254   
0255   m_error = beta/normRhs; 
0256   if(m_error < m_tolerance)
0257     m_info = Success; 
0258   else
0259     m_info = NoConvergence;
0260   
0261   // Iterative process
0262   while (nbIts < m_iterations && m_info == NoConvergence)
0263   {
0264     dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts); 
0265     
0266     // Compute the new residual vector for the restart 
0267     if (nbIts < m_iterations && m_info == NoConvergence) {
0268       r0 = rhs - mat * x;
0269       beta = r0.norm();
0270     }
0271   }
0272 } 
0273 
0274 /**
0275  * \brief Perform one restart cycle of DGMRES
0276  * \param mat The coefficient matrix
0277  * \param precond The preconditioner
0278  * \param x the new approximated solution
0279  * \param r0 The initial residual vector
0280  * \param beta The norm of the residual computed so far
0281  * \param normRhs The norm of the right hand side vector
0282  * \param nbIts The number of iterations
0283  */
0284 template< typename _MatrixType, typename _Preconditioner>
0285 template<typename Dest>
0286 Index DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, Index& nbIts) const
0287 {
0288   //Initialization 
0289   DenseVector g(m_restart+1); // Right hand side of the least square problem
0290   g.setZero();  
0291   g(0) = Scalar(beta); 
0292   m_V.col(0) = r0/beta; 
0293   m_info = NoConvergence; 
0294   std::vector<JacobiRotation<Scalar> >gr(m_restart); // Givens rotations
0295   Index it = 0; // Number of inner iterations 
0296   Index n = mat.rows();
0297   DenseVector tv1(n), tv2(n);  //Temporary vectors
0298   while (m_info == NoConvergence && it < m_restart && nbIts < m_iterations)
0299   {    
0300     // Apply preconditioner(s) at right
0301     if (m_isDeflInitialized )
0302     {
0303       dgmresApplyDeflation(m_V.col(it), tv1); // Deflation
0304       tv2 = precond.solve(tv1); 
0305     }
0306     else
0307     {
0308       tv2 = precond.solve(m_V.col(it)); // User's selected preconditioner
0309     }
0310     tv1 = mat * tv2; 
0311    
0312     // Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt
0313     Scalar coef; 
0314     for (Index i = 0; i <= it; ++i)
0315     { 
0316       coef = tv1.dot(m_V.col(i));
0317       tv1 = tv1 - coef * m_V.col(i); 
0318       m_H(i,it) = coef; 
0319       m_Hes(i,it) = coef; 
0320     }
0321     // Normalize the vector 
0322     coef = tv1.norm(); 
0323     m_V.col(it+1) = tv1/coef;
0324     m_H(it+1, it) = coef;
0325 //     m_Hes(it+1,it) = coef; 
0326     
0327     // FIXME Check for happy breakdown 
0328     
0329     // Update Hessenberg matrix with Givens rotations
0330     for (Index i = 1; i <= it; ++i) 
0331     {
0332       m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint());
0333     }
0334     // Compute the new plane rotation 
0335     gr[it].makeGivens(m_H(it, it), m_H(it+1,it)); 
0336     // Apply the new rotation
0337     m_H.col(it).applyOnTheLeft(it,it+1,gr[it].adjoint());
0338     g.applyOnTheLeft(it,it+1, gr[it].adjoint()); 
0339     
0340     beta = std::abs(g(it+1));
0341     m_error = beta/normRhs; 
0342     // std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl;
0343     it++; nbIts++; 
0344     
0345     if (m_error < m_tolerance)
0346     {
0347       // The method has converged
0348       m_info = Success;
0349       break;
0350     }
0351   }
0352   
0353   // Compute the new coefficients by solving the least square problem
0354 //   it++;
0355   //FIXME  Check first if the matrix is singular ... zero diagonal
0356   DenseVector nrs(m_restart); 
0357   nrs = m_H.topLeftCorner(it,it).template triangularView<Upper>().solve(g.head(it)); 
0358   
0359   // Form the new solution
0360   if (m_isDeflInitialized)
0361   {
0362     tv1 = m_V.leftCols(it) * nrs; 
0363     dgmresApplyDeflation(tv1, tv2); 
0364     x = x + precond.solve(tv2);
0365   }
0366   else
0367     x = x + precond.solve(m_V.leftCols(it) * nrs); 
0368   
0369   // Go for a new cycle and compute data for deflation
0370   if(nbIts < m_iterations && m_info == NoConvergence && m_neig > 0 && (m_r+m_neig) < m_maxNeig)
0371     dgmresComputeDeflationData(mat, precond, it, m_neig); 
0372   return 0; 
0373   
0374 }
0375 
0376 
0377 template< typename _MatrixType, typename _Preconditioner>
0378 void DGMRES<_MatrixType, _Preconditioner>::dgmresInitDeflation(Index& rows) const
0379 {
0380   m_U.resize(rows, m_maxNeig);
0381   m_MU.resize(rows, m_maxNeig); 
0382   m_T.resize(m_maxNeig, m_maxNeig);
0383   m_lambdaN = 0.0; 
0384   m_isDeflAllocated = true; 
0385 }
0386 
0387 template< typename _MatrixType, typename _Preconditioner>
0388 inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const ComplexSchur<DenseMatrix>& schurofH) const
0389 {
0390   return schurofH.matrixT().diagonal();
0391 }
0392 
0393 template< typename _MatrixType, typename _Preconditioner>
0394 inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const RealSchur<DenseMatrix>& schurofH) const
0395 {
0396   const DenseMatrix& T = schurofH.matrixT();
0397   Index it = T.rows();
0398   ComplexVector eig(it);
0399   Index j = 0;
0400   while (j < it-1)
0401   {
0402     if (T(j+1,j) ==Scalar(0))
0403     {
0404       eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0)); 
0405       j++; 
0406     }
0407     else
0408     {
0409       eig(j) = std::complex<RealScalar>(T(j,j),T(j+1,j)); 
0410       eig(j+1) = std::complex<RealScalar>(T(j,j+1),T(j+1,j+1));
0411       j++;
0412     }
0413   }
0414   if (j < it-1) eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0));
0415   return eig;
0416 }
0417 
0418 template< typename _MatrixType, typename _Preconditioner>
0419 Index DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const
0420 {
0421   // First, find the Schur form of the Hessenberg matrix H
0422   typename internal::conditional<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> >::type schurofH; 
0423   bool computeU = true;
0424   DenseMatrix matrixQ(it,it); 
0425   matrixQ.setIdentity();
0426   schurofH.computeFromHessenberg(m_Hes.topLeftCorner(it,it), matrixQ, computeU); 
0427   
0428   ComplexVector eig(it);
0429   Matrix<StorageIndex,Dynamic,1>perm(it);
0430   eig = this->schurValues(schurofH);
0431   
0432   // Reorder the absolute values of Schur values
0433   DenseRealVector modulEig(it); 
0434   for (Index j=0; j<it; ++j) modulEig(j) = std::abs(eig(j)); 
0435   perm.setLinSpaced(it,0,internal::convert_index<StorageIndex>(it-1));
0436   internal::sortWithPermutation(modulEig, perm, neig);
0437   
0438   if (!m_lambdaN)
0439   {
0440     m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN);
0441   }
0442   //Count the real number of extracted eigenvalues (with complex conjugates)
0443   Index nbrEig = 0; 
0444   while (nbrEig < neig)
0445   {
0446     if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++; 
0447     else nbrEig += 2; 
0448   }
0449   // Extract the  Schur vectors corresponding to the smallest Ritz values
0450   DenseMatrix Sr(it, nbrEig); 
0451   Sr.setZero();
0452   for (Index j = 0; j < nbrEig; j++)
0453   {
0454     Sr.col(j) = schurofH.matrixU().col(perm(it-j-1));
0455   }
0456   
0457   // Form the Schur vectors of the initial matrix using the Krylov basis
0458   DenseMatrix X; 
0459   X = m_V.leftCols(it) * Sr;
0460   if (m_r)
0461   {
0462    // Orthogonalize X against m_U using modified Gram-Schmidt
0463    for (Index j = 0; j < nbrEig; j++)
0464      for (Index k =0; k < m_r; k++)
0465       X.col(j) = X.col(j) - (m_U.col(k).dot(X.col(j)))*m_U.col(k); 
0466   }
0467   
0468   // Compute m_MX = A * M^-1 * X
0469   Index m = m_V.rows();
0470   if (!m_isDeflAllocated) 
0471     dgmresInitDeflation(m); 
0472   DenseMatrix MX(m, nbrEig);
0473   DenseVector tv1(m);
0474   for (Index j = 0; j < nbrEig; j++)
0475   {
0476     tv1 = mat * X.col(j);
0477     MX.col(j) = precond.solve(tv1);
0478   }
0479   
0480   //Update m_T = [U'MU U'MX; X'MU X'MX]
0481   m_T.block(m_r, m_r, nbrEig, nbrEig) = X.transpose() * MX; 
0482   if(m_r)
0483   {
0484     m_T.block(0, m_r, m_r, nbrEig) = m_U.leftCols(m_r).transpose() * MX; 
0485     m_T.block(m_r, 0, nbrEig, m_r) = X.transpose() * m_MU.leftCols(m_r);
0486   }
0487   
0488   // Save X into m_U and m_MX in m_MU
0489   for (Index j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j);
0490   for (Index j = 0; j < nbrEig; j++) m_MU.col(m_r+j) = MX.col(j);
0491   // Increase the size of the invariant subspace
0492   m_r += nbrEig; 
0493   
0494   // Factorize m_T into m_luT
0495   m_luT.compute(m_T.topLeftCorner(m_r, m_r));
0496   
0497   //FIXME CHeck if the factorization was correctly done (nonsingular matrix)
0498   m_isDeflInitialized = true;
0499   return 0; 
0500 }
0501 template<typename _MatrixType, typename _Preconditioner>
0502 template<typename RhsType, typename DestType>
0503 Index DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x, DestType &y) const
0504 {
0505   DenseVector x1 = m_U.leftCols(m_r).transpose() * x; 
0506   y = x + m_U.leftCols(m_r) * ( m_lambdaN * m_luT.solve(x1) - x1);
0507   return 0; 
0508 }
0509 
0510 } // end namespace Eigen
0511 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!