$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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