![]() |
Eigen-unsupported
3.3.3
|
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 bool flag; 00044 for (Index k = 0; k < ncut; k++) 00045 { 00046 flag = false; 00047 for (Index j = 0; j < vec.size()-1; j++) 00048 { 00049 if ( vec(perm(j)) < vec(perm(j+1)) ) 00050 { 00051 std::swap(perm(j),perm(j+1)); 00052 flag = true; 00053 } 00054 if (!flag) break; // The vector is in sorted order 00055 } 00056 } 00057 } 00058 00059 } 00101 template< typename _MatrixType, typename _Preconditioner> 00102 class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > 00103 { 00104 typedef IterativeSolverBase<DGMRES> Base; 00105 using Base::matrix; 00106 using Base::m_error; 00107 using Base::m_iterations; 00108 using Base::m_info; 00109 using Base::m_isInitialized; 00110 using Base::m_tolerance; 00111 public: 00112 using Base::_solve_impl; 00113 typedef _MatrixType MatrixType; 00114 typedef typename MatrixType::Scalar Scalar; 00115 typedef typename MatrixType::Index Index; 00116 typedef typename MatrixType::StorageIndex StorageIndex; 00117 typedef typename MatrixType::RealScalar RealScalar; 00118 typedef _Preconditioner Preconditioner; 00119 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 00120 typedef Matrix<RealScalar,Dynamic,Dynamic> DenseRealMatrix; 00121 typedef Matrix<Scalar,Dynamic,1> DenseVector; 00122 typedef Matrix<RealScalar,Dynamic,1> DenseRealVector; 00123 typedef Matrix<std::complex<RealScalar>, Dynamic, 1> ComplexVector; 00124 00125 00127 DGMRES() : Base(),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {} 00128 00139 template<typename MatrixDerived> 00140 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) {} 00141 00142 ~DGMRES() {} 00143 00145 template<typename Rhs,typename Dest> 00146 void _solve_with_guess_impl(const Rhs& b, Dest& x) const 00147 { 00148 bool failed = false; 00149 for(int j=0; j<b.cols(); ++j) 00150 { 00151 m_iterations = Base::maxIterations(); 00152 m_error = Base::m_tolerance; 00153 00154 typename Dest::ColXpr xj(x,j); 00155 dgmres(matrix(), b.col(j), xj, Base::m_preconditioner); 00156 } 00157 m_info = failed ? NumericalIssue 00158 : m_error <= Base::m_tolerance ? Success 00159 : NoConvergence; 00160 m_isInitialized = true; 00161 } 00162 00164 template<typename Rhs,typename Dest> 00165 void _solve_impl(const Rhs& b, MatrixBase<Dest>& x) const 00166 { 00167 x = b; 00168 _solve_with_guess_impl(b,x.derived()); 00169 } 00173 int restart() { return m_restart; } 00174 00178 void set_restart(const int restart) { m_restart=restart; } 00179 00183 void setEigenv(const int neig) 00184 { 00185 m_neig = neig; 00186 if (neig+1 > m_maxNeig) m_maxNeig = neig+1; // To allow for complex conjugates 00187 } 00188 00192 int deflSize() {return m_r; } 00193 00197 void setMaxEigenv(const int maxNeig) { m_maxNeig = maxNeig; } 00198 00199 protected: 00200 // DGMRES algorithm 00201 template<typename Rhs, typename Dest> 00202 void dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, const Preconditioner& precond) const; 00203 // Perform one cycle of GMRES 00204 template<typename Dest> 00205 int dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const; 00206 // Compute data to use for deflation 00207 int dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const; 00208 // Apply deflation to a vector 00209 template<typename RhsType, typename DestType> 00210 int dgmresApplyDeflation(const RhsType& In, DestType& Out) const; 00211 ComplexVector schurValues(const ComplexSchur<DenseMatrix>& schurofH) const; 00212 ComplexVector schurValues(const RealSchur<DenseMatrix>& schurofH) const; 00213 // Init data for deflation 00214 void dgmresInitDeflation(Index& rows) const; 00215 mutable DenseMatrix m_V; // Krylov basis vectors 00216 mutable DenseMatrix m_H; // Hessenberg matrix 00217 mutable DenseMatrix m_Hes; // Initial hessenberg matrix wihout Givens rotations applied 00218 mutable Index m_restart; // Maximum size of the Krylov subspace 00219 mutable DenseMatrix m_U; // Vectors that form the basis of the invariant subspace 00220 mutable DenseMatrix m_MU; // matrix operator applied to m_U (for next cycles) 00221 mutable DenseMatrix m_T; /* T=U^T*M^{-1}*A*U */ 00222 mutable PartialPivLU<DenseMatrix> m_luT; // LU factorization of m_T 00223 mutable StorageIndex m_neig; //Number of eigenvalues to extract at each restart 00224 mutable int m_r; // Current number of deflated eigenvalues, size of m_U 00225 mutable int m_maxNeig; // Maximum number of eigenvalues to deflate 00226 mutable RealScalar m_lambdaN; //Modulus of the largest eigenvalue of A 00227 mutable bool m_isDeflAllocated; 00228 mutable bool m_isDeflInitialized; 00229 00230 //Adaptive strategy 00231 mutable RealScalar m_smv; // Smaller multiple of the remaining number of steps allowed 00232 mutable bool m_force; // Force the use of deflation at each restart 00233 00234 }; 00241 template< typename _MatrixType, typename _Preconditioner> 00242 template<typename Rhs, typename Dest> 00243 void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, 00244 const Preconditioner& precond) const 00245 { 00246 //Initialization 00247 int n = mat.rows(); 00248 DenseVector r0(n); 00249 int nbIts = 0; 00250 m_H.resize(m_restart+1, m_restart); 00251 m_Hes.resize(m_restart, m_restart); 00252 m_V.resize(n,m_restart+1); 00253 //Initial residual vector and intial norm 00254 x = precond.solve(x); 00255 r0 = rhs - mat * x; 00256 RealScalar beta = r0.norm(); 00257 RealScalar normRhs = rhs.norm(); 00258 m_error = beta/normRhs; 00259 if(m_error < m_tolerance) 00260 m_info = Success; 00261 else 00262 m_info = NoConvergence; 00263 00264 // Iterative process 00265 while (nbIts < m_iterations && m_info == NoConvergence) 00266 { 00267 dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts); 00268 00269 // Compute the new residual vector for the restart 00270 if (nbIts < m_iterations && m_info == NoConvergence) 00271 r0 = rhs - mat * x; 00272 } 00273 } 00274 00285 template< typename _MatrixType, typename _Preconditioner> 00286 template<typename Dest> 00287 int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const 00288 { 00289 //Initialization 00290 DenseVector g(m_restart+1); // Right hand side of the least square problem 00291 g.setZero(); 00292 g(0) = Scalar(beta); 00293 m_V.col(0) = r0/beta; 00294 m_info = NoConvergence; 00295 std::vector<JacobiRotation<Scalar> >gr(m_restart); // Givens rotations 00296 int it = 0; // Number of inner iterations 00297 int n = mat.rows(); 00298 DenseVector tv1(n), tv2(n); //Temporary vectors 00299 while (m_info == NoConvergence && it < m_restart && nbIts < m_iterations) 00300 { 00301 // Apply preconditioner(s) at right 00302 if (m_isDeflInitialized ) 00303 { 00304 dgmresApplyDeflation(m_V.col(it), tv1); // Deflation 00305 tv2 = precond.solve(tv1); 00306 } 00307 else 00308 { 00309 tv2 = precond.solve(m_V.col(it)); // User's selected preconditioner 00310 } 00311 tv1 = mat * tv2; 00312 00313 // Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt 00314 Scalar coef; 00315 for (int i = 0; i <= it; ++i) 00316 { 00317 coef = tv1.dot(m_V.col(i)); 00318 tv1 = tv1 - coef * m_V.col(i); 00319 m_H(i,it) = coef; 00320 m_Hes(i,it) = coef; 00321 } 00322 // Normalize the vector 00323 coef = tv1.norm(); 00324 m_V.col(it+1) = tv1/coef; 00325 m_H(it+1, it) = coef; 00326 // m_Hes(it+1,it) = coef; 00327 00328 // FIXME Check for happy breakdown 00329 00330 // Update Hessenberg matrix with Givens rotations 00331 for (int i = 1; i <= it; ++i) 00332 { 00333 m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint()); 00334 } 00335 // Compute the new plane rotation 00336 gr[it].makeGivens(m_H(it, it), m_H(it+1,it)); 00337 // Apply the new rotation 00338 m_H.col(it).applyOnTheLeft(it,it+1,gr[it].adjoint()); 00339 g.applyOnTheLeft(it,it+1, gr[it].adjoint()); 00340 00341 beta = std::abs(g(it+1)); 00342 m_error = beta/normRhs; 00343 // std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl; 00344 it++; nbIts++; 00345 00346 if (m_error < m_tolerance) 00347 { 00348 // The method has converged 00349 m_info = Success; 00350 break; 00351 } 00352 } 00353 00354 // Compute the new coefficients by solving the least square problem 00355 // it++; 00356 //FIXME Check first if the matrix is singular ... zero diagonal 00357 DenseVector nrs(m_restart); 00358 nrs = m_H.topLeftCorner(it,it).template triangularView<Upper>().solve(g.head(it)); 00359 00360 // Form the new solution 00361 if (m_isDeflInitialized) 00362 { 00363 tv1 = m_V.leftCols(it) * nrs; 00364 dgmresApplyDeflation(tv1, tv2); 00365 x = x + precond.solve(tv2); 00366 } 00367 else 00368 x = x + precond.solve(m_V.leftCols(it) * nrs); 00369 00370 // Go for a new cycle and compute data for deflation 00371 if(nbIts < m_iterations && m_info == NoConvergence && m_neig > 0 && (m_r+m_neig) < m_maxNeig) 00372 dgmresComputeDeflationData(mat, precond, it, m_neig); 00373 return 0; 00374 00375 } 00376 00377 00378 template< typename _MatrixType, typename _Preconditioner> 00379 void DGMRES<_MatrixType, _Preconditioner>::dgmresInitDeflation(Index& rows) const 00380 { 00381 m_U.resize(rows, m_maxNeig); 00382 m_MU.resize(rows, m_maxNeig); 00383 m_T.resize(m_maxNeig, m_maxNeig); 00384 m_lambdaN = 0.0; 00385 m_isDeflAllocated = true; 00386 } 00387 00388 template< typename _MatrixType, typename _Preconditioner> 00389 inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const ComplexSchur<DenseMatrix>& schurofH) const 00390 { 00391 return schurofH.matrixT().diagonal(); 00392 } 00393 00394 template< typename _MatrixType, typename _Preconditioner> 00395 inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const RealSchur<DenseMatrix>& schurofH) const 00396 { 00397 typedef typename MatrixType::Index Index; 00398 const DenseMatrix& T = schurofH.matrixT(); 00399 Index it = T.rows(); 00400 ComplexVector eig(it); 00401 Index j = 0; 00402 while (j < it-1) 00403 { 00404 if (T(j+1,j) ==Scalar(0)) 00405 { 00406 eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0)); 00407 j++; 00408 } 00409 else 00410 { 00411 eig(j) = std::complex<RealScalar>(T(j,j),T(j+1,j)); 00412 eig(j+1) = std::complex<RealScalar>(T(j,j+1),T(j+1,j+1)); 00413 j++; 00414 } 00415 } 00416 if (j < it-1) eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0)); 00417 return eig; 00418 } 00419 00420 template< typename _MatrixType, typename _Preconditioner> 00421 int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const 00422 { 00423 // First, find the Schur form of the Hessenberg matrix H 00424 typename internal::conditional<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> >::type schurofH; 00425 bool computeU = true; 00426 DenseMatrix matrixQ(it,it); 00427 matrixQ.setIdentity(); 00428 schurofH.computeFromHessenberg(m_Hes.topLeftCorner(it,it), matrixQ, computeU); 00429 00430 ComplexVector eig(it); 00431 Matrix<StorageIndex,Dynamic,1>perm(it); 00432 eig = this->schurValues(schurofH); 00433 00434 // Reorder the absolute values of Schur values 00435 DenseRealVector modulEig(it); 00436 for (int j=0; j<it; ++j) modulEig(j) = std::abs(eig(j)); 00437 perm.setLinSpaced(it,0,it-1); 00438 internal::sortWithPermutation(modulEig, perm, neig); 00439 00440 if (!m_lambdaN) 00441 { 00442 m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN); 00443 } 00444 //Count the real number of extracted eigenvalues (with complex conjugates) 00445 int nbrEig = 0; 00446 while (nbrEig < neig) 00447 { 00448 if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++; 00449 else nbrEig += 2; 00450 } 00451 // Extract the Schur vectors corresponding to the smallest Ritz values 00452 DenseMatrix Sr(it, nbrEig); 00453 Sr.setZero(); 00454 for (int j = 0; j < nbrEig; j++) 00455 { 00456 Sr.col(j) = schurofH.matrixU().col(perm(it-j-1)); 00457 } 00458 00459 // Form the Schur vectors of the initial matrix using the Krylov basis 00460 DenseMatrix X; 00461 X = m_V.leftCols(it) * Sr; 00462 if (m_r) 00463 { 00464 // Orthogonalize X against m_U using modified Gram-Schmidt 00465 for (int j = 0; j < nbrEig; j++) 00466 for (int k =0; k < m_r; k++) 00467 X.col(j) = X.col(j) - (m_U.col(k).dot(X.col(j)))*m_U.col(k); 00468 } 00469 00470 // Compute m_MX = A * M^-1 * X 00471 Index m = m_V.rows(); 00472 if (!m_isDeflAllocated) 00473 dgmresInitDeflation(m); 00474 DenseMatrix MX(m, nbrEig); 00475 DenseVector tv1(m); 00476 for (int j = 0; j < nbrEig; j++) 00477 { 00478 tv1 = mat * X.col(j); 00479 MX.col(j) = precond.solve(tv1); 00480 } 00481 00482 //Update m_T = [U'MU U'MX; X'MU X'MX] 00483 m_T.block(m_r, m_r, nbrEig, nbrEig) = X.transpose() * MX; 00484 if(m_r) 00485 { 00486 m_T.block(0, m_r, m_r, nbrEig) = m_U.leftCols(m_r).transpose() * MX; 00487 m_T.block(m_r, 0, nbrEig, m_r) = X.transpose() * m_MU.leftCols(m_r); 00488 } 00489 00490 // Save X into m_U and m_MX in m_MU 00491 for (int j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j); 00492 for (int j = 0; j < nbrEig; j++) m_MU.col(m_r+j) = MX.col(j); 00493 // Increase the size of the invariant subspace 00494 m_r += nbrEig; 00495 00496 // Factorize m_T into m_luT 00497 m_luT.compute(m_T.topLeftCorner(m_r, m_r)); 00498 00499 //FIXME CHeck if the factorization was correctly done (nonsingular matrix) 00500 m_isDeflInitialized = true; 00501 return 0; 00502 } 00503 template<typename _MatrixType, typename _Preconditioner> 00504 template<typename RhsType, typename DestType> 00505 int DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x, DestType &y) const 00506 { 00507 DenseVector x1 = m_U.leftCols(m_r).transpose() * x; 00508 y = x + m_U.leftCols(m_r) * ( m_lambdaN * m_luT.solve(x1) - x1); 00509 return 0; 00510 } 00511 00512 } // end namespace Eigen 00513 #endif