DGMRES.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_DGMRES_H
00011 #define EIGEN_DGMRES_H
00012 
00013 #include <Eigen/Eigenvalues>
00014 
00015 namespace Eigen { 
00016   
00017 template< typename _MatrixType,
00018           typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
00019 class DGMRES;
00020 
00021 namespace internal {
00022 
00023 template< typename _MatrixType, typename _Preconditioner>
00024 struct traits<DGMRES<_MatrixType,_Preconditioner> >
00025 {
00026   typedef _MatrixType MatrixType;
00027   typedef _Preconditioner Preconditioner;
00028 };
00029 
00038 template <typename VectorType, typename IndexType>
00039 void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType::Scalar& ncut)
00040 {
00041   eigen_assert(vec.size() == perm.size());
00042   typedef typename IndexType::Scalar Index; 
00043   typedef typename VectorType::Scalar Scalar; 
00044   bool flag; 
00045   for (Index k  = 0; k < ncut; k++)
00046   {
00047     flag = false;
00048     for (Index j = 0; j < vec.size()-1; j++)
00049     {
00050       if ( vec(perm(j)) < vec(perm(j+1)) )
00051       {
00052         std::swap(perm(j),perm(j+1)); 
00053         flag = true;
00054       }
00055       if (!flag) break; // The vector is in sorted order
00056     }
00057   }
00058 }
00059 
00060 }
00100 template< typename _MatrixType, typename _Preconditioner>
00101 class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
00102 {
00103     typedef IterativeSolverBase<DGMRES> Base;
00104     using Base::mp_matrix;
00105     using Base::m_error;
00106     using Base::m_iterations;
00107     using Base::m_info;
00108     using Base::m_isInitialized;
00109     using Base::m_tolerance; 
00110   public:
00111     typedef _MatrixType MatrixType;
00112     typedef typename MatrixType::Scalar Scalar;
00113     typedef typename MatrixType::Index Index;
00114     typedef typename MatrixType::RealScalar RealScalar;
00115     typedef _Preconditioner Preconditioner;
00116     typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 
00117     typedef Matrix<RealScalar,Dynamic,Dynamic> DenseRealMatrix; 
00118     typedef Matrix<Scalar,Dynamic,1> DenseVector;
00119     typedef Matrix<RealScalar,Dynamic,1> DenseRealVector; 
00120     typedef Matrix<std::complex<RealScalar>, Dynamic, 1> ComplexVector;
00121  
00122     
00124   DGMRES() : Base(),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {}
00125 
00136   DGMRES(const MatrixType& A) : Base(A),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false)
00137   {}
00138 
00139   ~DGMRES() {}
00140   
00146   template<typename Rhs,typename Guess>
00147   inline const internal::solve_retval_with_guess<DGMRES, Rhs, Guess>
00148   solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
00149   {
00150     eigen_assert(m_isInitialized && "DGMRES is not initialized.");
00151     eigen_assert(Base::rows()==b.rows()
00152               && "DGMRES::solve(): invalid number of rows of the right hand side matrix b");
00153     return internal::solve_retval_with_guess
00154             <DGMRES, Rhs, Guess>(*this, b.derived(), x0);
00155   }
00156   
00158   template<typename Rhs,typename Dest>
00159   void _solveWithGuess(const Rhs& b, Dest& x) const
00160   {    
00161     bool failed = false;
00162     for(int j=0; j<b.cols(); ++j)
00163     {
00164       m_iterations = Base::maxIterations();
00165       m_error = Base::m_tolerance;
00166       
00167       typename Dest::ColXpr xj(x,j);
00168       dgmres(*mp_matrix, b.col(j), xj, Base::m_preconditioner);
00169     }
00170     m_info = failed ? NumericalIssue
00171            : m_error <= Base::m_tolerance ? Success
00172            : NoConvergence;
00173     m_isInitialized = true;
00174   }
00175 
00177   template<typename Rhs,typename Dest>
00178   void _solve(const Rhs& b, Dest& x) const
00179   {
00180     x = b;
00181     _solveWithGuess(b,x);
00182   }
00186   int restart() { return m_restart; }
00187   
00191   void set_restart(const int restart) { m_restart=restart; }
00192   
00196   void setEigenv(const int neig) 
00197   {
00198     m_neig = neig;
00199     if (neig+1 > m_maxNeig) m_maxNeig = neig+1; // To allow for complex conjugates
00200   }
00201   
00205   int deflSize() {return m_r; }
00206   
00210   void setMaxEigenv(const int maxNeig) { m_maxNeig = maxNeig; }
00211   
00212   protected:
00213     // DGMRES algorithm 
00214     template<typename Rhs, typename Dest>
00215     void dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, const Preconditioner& precond) const;
00216     // Perform one cycle of GMRES
00217     template<typename Dest>
00218     int dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const; 
00219     // Compute data to use for deflation 
00220     int dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, Index& neig) const;
00221     // Apply deflation to a vector
00222     template<typename RhsType, typename DestType>
00223     int dgmresApplyDeflation(const RhsType& In, DestType& Out) const; 
00224     ComplexVector schurValues(const ComplexSchur<DenseMatrix>& schurofH) const;
00225     ComplexVector schurValues(const RealSchur<DenseMatrix>& schurofH) const;
00226     // Init data for deflation
00227     void dgmresInitDeflation(Index& rows) const; 
00228     mutable DenseMatrix m_V; // Krylov basis vectors
00229     mutable DenseMatrix m_H; // Hessenberg matrix 
00230     mutable DenseMatrix m_Hes; // Initial hessenberg matrix wihout Givens rotations applied
00231     mutable Index m_restart; // Maximum size of the Krylov subspace
00232     mutable DenseMatrix m_U; // Vectors that form the basis of the invariant subspace 
00233     mutable DenseMatrix m_MU; // matrix operator applied to m_U (for next cycles)
00234     mutable DenseMatrix m_T; /* T=U^T*M^{-1}*A*U */
00235     mutable PartialPivLU<DenseMatrix> m_luT; // LU factorization of m_T
00236     mutable int m_neig; //Number of eigenvalues to extract at each restart
00237     mutable int m_r; // Current number of deflated eigenvalues, size of m_U
00238     mutable int m_maxNeig; // Maximum number of eigenvalues to deflate
00239     mutable RealScalar m_lambdaN; //Modulus of the largest eigenvalue of A
00240     mutable bool m_isDeflAllocated;
00241     mutable bool m_isDeflInitialized;
00242     
00243     //Adaptive strategy 
00244     mutable RealScalar m_smv; // Smaller multiple of the remaining number of steps allowed
00245     mutable bool m_force; // Force the use of deflation at each restart
00246     
00247 }; 
00254 template< typename _MatrixType, typename _Preconditioner>
00255 template<typename Rhs, typename Dest>
00256 void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x,
00257               const Preconditioner& precond) const
00258 {
00259   //Initialization
00260   int n = mat.rows(); 
00261   DenseVector r0(n); 
00262   int nbIts = 0; 
00263   m_H.resize(m_restart+1, m_restart);
00264   m_Hes.resize(m_restart, m_restart);
00265   m_V.resize(n,m_restart+1);
00266   //Initial residual vector and intial norm
00267   x = precond.solve(x);
00268   r0 = rhs - mat * x; 
00269   RealScalar beta = r0.norm(); 
00270   RealScalar normRhs = rhs.norm();
00271   m_error = beta/normRhs; 
00272   if(m_error < m_tolerance)
00273     m_info = Success; 
00274   else
00275     m_info = NoConvergence;
00276   
00277   // Iterative process
00278   while (nbIts < m_iterations && m_info == NoConvergence)
00279   {
00280     dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts); 
00281     
00282     // Compute the new residual vector for the restart 
00283     if (nbIts < m_iterations && m_info == NoConvergence)
00284       r0 = rhs - mat * x; 
00285   }
00286 } 
00287 
00298 template< typename _MatrixType, typename _Preconditioner>
00299 template<typename Dest>
00300 int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const
00301 {
00302   //Initialization 
00303   DenseVector g(m_restart+1); // Right hand side of the least square problem
00304   g.setZero();  
00305   g(0) = Scalar(beta); 
00306   m_V.col(0) = r0/beta; 
00307   m_info = NoConvergence; 
00308   std::vector<JacobiRotation<Scalar> >gr(m_restart); // Givens rotations
00309   int it = 0; // Number of inner iterations 
00310   int n = mat.rows();
00311   DenseVector tv1(n), tv2(n);  //Temporary vectors
00312   while (m_info == NoConvergence && it < m_restart && nbIts < m_iterations)
00313   {    
00314     // Apply preconditioner(s) at right
00315     if (m_isDeflInitialized )
00316     {
00317       dgmresApplyDeflation(m_V.col(it), tv1); // Deflation
00318       tv2 = precond.solve(tv1); 
00319     }
00320     else
00321     {
00322       tv2 = precond.solve(m_V.col(it)); // User's selected preconditioner
00323     }
00324     tv1 = mat * tv2; 
00325    
00326     // Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt
00327     Scalar coef; 
00328     for (int i = 0; i <= it; ++i)
00329     { 
00330       coef = tv1.dot(m_V.col(i));
00331       tv1 = tv1 - coef * m_V.col(i); 
00332       m_H(i,it) = coef; 
00333       m_Hes(i,it) = coef; 
00334     }
00335     // Normalize the vector 
00336     coef = tv1.norm(); 
00337     m_V.col(it+1) = tv1/coef;
00338     m_H(it+1, it) = coef;
00339 //     m_Hes(it+1,it) = coef; 
00340     
00341     // FIXME Check for happy breakdown 
00342     
00343     // Update Hessenberg matrix with Givens rotations
00344     for (int i = 1; i <= it; ++i) 
00345     {
00346       m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint());
00347     }
00348     // Compute the new plane rotation 
00349     gr[it].makeGivens(m_H(it, it), m_H(it+1,it)); 
00350     // Apply the new rotation
00351     m_H.col(it).applyOnTheLeft(it,it+1,gr[it].adjoint());
00352     g.applyOnTheLeft(it,it+1, gr[it].adjoint()); 
00353     
00354     beta = std::abs(g(it+1));
00355     m_error = beta/normRhs; 
00356     std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl;
00357     it++; nbIts++; 
00358     
00359     if (m_error < m_tolerance)
00360     {
00361       // The method has converged
00362       m_info = Success;
00363       break;
00364     }
00365   }
00366   
00367   // Compute the new coefficients by solving the least square problem
00368 //   it++;
00369   //FIXME  Check first if the matrix is singular ... zero diagonal
00370   DenseVector nrs(m_restart); 
00371   nrs = m_H.topLeftCorner(it,it).template triangularView<Upper>().solve(g.head(it)); 
00372   
00373   // Form the new solution
00374   if (m_isDeflInitialized)
00375   {
00376     tv1 = m_V.leftCols(it) * nrs; 
00377     dgmresApplyDeflation(tv1, tv2); 
00378     x = x + precond.solve(tv2);
00379   }
00380   else
00381     x = x + precond.solve(m_V.leftCols(it) * nrs); 
00382   
00383   // Go for a new cycle and compute data for deflation
00384   if(nbIts < m_iterations && m_info == NoConvergence && m_neig > 0 && (m_r+m_neig) < m_maxNeig)
00385     dgmresComputeDeflationData(mat, precond, it, m_neig); 
00386   return 0; 
00387   
00388 }
00389 
00390 
00391 template< typename _MatrixType, typename _Preconditioner>
00392 void DGMRES<_MatrixType, _Preconditioner>::dgmresInitDeflation(Index& rows) const
00393 {
00394   m_U.resize(rows, m_maxNeig);
00395   m_MU.resize(rows, m_maxNeig); 
00396   m_T.resize(m_maxNeig, m_maxNeig);
00397   m_lambdaN = 0.0; 
00398   m_isDeflAllocated = true; 
00399 }
00400 
00401 template< typename _MatrixType, typename _Preconditioner>
00402 inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const ComplexSchur<DenseMatrix>& schurofH) const
00403 {
00404   return schurofH.matrixT().diagonal();
00405 }
00406 
00407 template< typename _MatrixType, typename _Preconditioner>
00408 inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const RealSchur<DenseMatrix>& schurofH) const
00409 {
00410   typedef typename MatrixType::Index Index;
00411   const DenseMatrix& T = schurofH.matrixT();
00412   Index it = T.rows();
00413   ComplexVector eig(it);
00414   Index j = 0;
00415   while (j < it-1)
00416   {
00417     if (T(j+1,j) ==Scalar(0))
00418     {
00419       eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0)); 
00420       j++; 
00421     }
00422     else
00423     {
00424       eig(j) = std::complex<RealScalar>(T(j,j),T(j+1,j)); 
00425       eig(j+1) = std::complex<RealScalar>(T(j,j+1),T(j+1,j+1));
00426       j++;
00427     }
00428   }
00429   if (j < it-1) eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0));
00430   return eig;
00431 }
00432 
00433 template< typename _MatrixType, typename _Preconditioner>
00434 int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, Index& neig) const
00435 {
00436   // First, find the Schur form of the Hessenberg matrix H
00437   typename internal::conditional<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> >::type schurofH; 
00438   bool computeU = true;
00439   DenseMatrix matrixQ(it,it); 
00440   matrixQ.setIdentity();
00441   schurofH.computeFromHessenberg(m_Hes.topLeftCorner(it,it), matrixQ, computeU); 
00442   
00443   ComplexVector eig(it);
00444   Matrix<Index,Dynamic,1>perm(it); 
00445   eig = this->schurValues(schurofH);
00446   
00447   // Reorder the absolute values of Schur values
00448   DenseRealVector modulEig(it); 
00449   for (int j=0; j<it; ++j) modulEig(j) = std::abs(eig(j)); 
00450   perm.setLinSpaced(it,0,it-1);
00451   internal::sortWithPermutation(modulEig, perm, neig);
00452   
00453   if (!m_lambdaN)
00454   {
00455     m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN);
00456   }
00457   //Count the real number of extracted eigenvalues (with complex conjugates)
00458   int nbrEig = 0; 
00459   while (nbrEig < neig)
00460   {
00461     if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++; 
00462     else nbrEig += 2; 
00463   }
00464   // Extract the  Schur vectors corresponding to the smallest Ritz values
00465   DenseMatrix Sr(it, nbrEig); 
00466   Sr.setZero();
00467   for (int j = 0; j < nbrEig; j++)
00468   {
00469     Sr.col(j) = schurofH.matrixU().col(perm(it-j-1));
00470   }
00471   
00472   // Form the Schur vectors of the initial matrix using the Krylov basis
00473   DenseMatrix X; 
00474   X = m_V.leftCols(it) * Sr;
00475   if (m_r)
00476   {
00477    // Orthogonalize X against m_U using modified Gram-Schmidt
00478    for (int j = 0; j < nbrEig; j++)
00479      for (int k =0; k < m_r; k++)
00480       X.col(j) = X.col(j) - (m_U.col(k).dot(X.col(j)))*m_U.col(k); 
00481   }
00482   
00483   // Compute m_MX = A * M^-1 * X
00484   Index m = m_V.rows();
00485   if (!m_isDeflAllocated) 
00486     dgmresInitDeflation(m); 
00487   DenseMatrix MX(m, nbrEig);
00488   DenseVector tv1(m);
00489   for (int j = 0; j < nbrEig; j++)
00490   {
00491     tv1 = mat * X.col(j);
00492     MX.col(j) = precond.solve(tv1);
00493   }
00494   
00495   //Update m_T = [U'MU U'MX; X'MU X'MX]
00496   m_T.block(m_r, m_r, nbrEig, nbrEig) = X.transpose() * MX; 
00497   if(m_r)
00498   {
00499     m_T.block(0, m_r, m_r, nbrEig) = m_U.leftCols(m_r).transpose() * MX; 
00500     m_T.block(m_r, 0, nbrEig, m_r) = X.transpose() * m_MU.leftCols(m_r);
00501   }
00502   
00503   // Save X into m_U and m_MX in m_MU
00504   for (int j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j);
00505   for (int j = 0; j < nbrEig; j++) m_MU.col(m_r+j) = MX.col(j);
00506   // Increase the size of the invariant subspace
00507   m_r += nbrEig; 
00508   
00509   // Factorize m_T into m_luT
00510   m_luT.compute(m_T.topLeftCorner(m_r, m_r));
00511   
00512   //FIXME CHeck if the factorization was correctly done (nonsingular matrix)
00513   m_isDeflInitialized = true;
00514   return 0; 
00515 }
00516 template<typename _MatrixType, typename _Preconditioner>
00517 template<typename RhsType, typename DestType>
00518 int DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x, DestType &y) const
00519 {
00520   DenseVector x1 = m_U.leftCols(m_r).transpose() * x; 
00521   y = x + m_U.leftCols(m_r) * ( m_lambdaN * m_luT.solve(x1) - x1);
00522   return 0; 
00523 }
00524 
00525 namespace internal {
00526 
00527   template<typename _MatrixType, typename _Preconditioner, typename Rhs>
00528 struct solve_retval<DGMRES<_MatrixType, _Preconditioner>, Rhs>
00529   : solve_retval_base<DGMRES<_MatrixType, _Preconditioner>, Rhs>
00530 {
00531   typedef DGMRES<_MatrixType, _Preconditioner> Dec;
00532   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00533 
00534   template<typename Dest> void evalTo(Dest& dst) const
00535   {
00536     dec()._solve(rhs(),dst);
00537   }
00538 };
00539 } // end namespace internal
00540 
00541 } // end namespace Eigen
00542 #endif 


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:36:58