MINRES.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 Giacomo Po <gpo@ucla.edu>
00005 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 
00012 #ifndef EIGEN_MINRES_H_
00013 #define EIGEN_MINRES_H_
00014 
00015 
00016 namespace Eigen {
00017     
00018     namespace internal {
00019         
00029         template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
00030         EIGEN_DONT_INLINE
00031         void minres(const MatrixType& mat, const Rhs& rhs, Dest& x,
00032                     const Preconditioner& precond, int& iters,
00033                     typename Dest::RealScalar& tol_error)
00034         {
00035             using std::sqrt;
00036             typedef typename Dest::RealScalar RealScalar;
00037             typedef typename Dest::Scalar Scalar;
00038             typedef Matrix<Scalar,Dynamic,1> VectorType;
00039 
00040             // initialize
00041             const int maxIters(iters);  // initialize maxIters to iters
00042             const int N(mat.cols());    // the size of the matrix
00043             const RealScalar rhsNorm2(rhs.squaredNorm());
00044             const RealScalar threshold2(tol_error*tol_error*rhsNorm2); // convergence threshold (compared to residualNorm2)
00045             
00046             // Initialize preconditioned Lanczos
00047 //            VectorType v_old(N); // will be initialized inside loop
00048             VectorType v( VectorType::Zero(N) ); //initialize v
00049             VectorType v_new(rhs-mat*x); //initialize v_new
00050             RealScalar residualNorm2(v_new.squaredNorm());
00051 //            VectorType w(N); // will be initialized inside loop
00052             VectorType w_new(precond.solve(v_new)); // initialize w_new
00053 //            RealScalar beta; // will be initialized inside loop
00054             RealScalar beta_new2(v_new.dot(w_new));
00055             eigen_assert(beta_new2 >= 0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
00056             RealScalar beta_new(sqrt(beta_new2));
00057             const RealScalar beta_one(beta_new);
00058             v_new /= beta_new;
00059             w_new /= beta_new;
00060             // Initialize other variables
00061             RealScalar c(1.0); // the cosine of the Givens rotation
00062             RealScalar c_old(1.0);
00063             RealScalar s(0.0); // the sine of the Givens rotation
00064             RealScalar s_old(0.0); // the sine of the Givens rotation
00065 //            VectorType p_oold(N); // will be initialized in loop
00066             VectorType p_old(VectorType::Zero(N)); // initialize p_old=0
00067             VectorType p(p_old); // initialize p=0
00068             RealScalar eta(1.0);
00069                         
00070             iters = 0; // reset iters
00071             while ( iters < maxIters ){
00072                 
00073                 // Preconditioned Lanczos
00074                 /* Note that there are 4 variants on the Lanczos algorithm. These are
00075                  * described in Paige, C. C. (1972). Computational variants of
00076                  * the Lanczos method for the eigenproblem. IMA Journal of Applied
00077                  * Mathematics, 10(3), 373–381. The current implementation corresponds 
00078                  * to the case A(2,7) in the paper. It also corresponds to 
00079                  * algorithm 6.14 in Y. Saad, Iterative Methods for Sparse Linear
00080                  * Systems, 2003 p.173. For the preconditioned version see 
00081                  * A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM (1987).
00082                  */
00083                 const RealScalar beta(beta_new);
00084 //                v_old = v; // update: at first time step, this makes v_old = 0 so value of beta doesn't matter
00085                 const VectorType v_old(v); // NOT SURE IF CREATING v_old EVERY ITERATION IS EFFICIENT
00086                 v = v_new; // update
00087 //                w = w_new; // update
00088                 const VectorType w(w_new); // NOT SURE IF CREATING w EVERY ITERATION IS EFFICIENT
00089                 v_new.noalias() = mat*w - beta*v_old; // compute v_new
00090                 const RealScalar alpha = v_new.dot(w);
00091                 v_new -= alpha*v; // overwrite v_new
00092                 w_new = precond.solve(v_new); // overwrite w_new
00093                 beta_new2 = v_new.dot(w_new); // compute beta_new
00094                 eigen_assert(beta_new2 >= 0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
00095                 beta_new = sqrt(beta_new2); // compute beta_new
00096                 v_new /= beta_new; // overwrite v_new for next iteration
00097                 w_new /= beta_new; // overwrite w_new for next iteration
00098                 
00099                 // Givens rotation
00100                 const RealScalar r2 =s*alpha+c*c_old*beta; // s, s_old, c and c_old are still from previous iteration
00101                 const RealScalar r3 =s_old*beta; // s, s_old, c and c_old are still from previous iteration
00102                 const RealScalar r1_hat=c*alpha-c_old*s*beta;
00103                 const RealScalar r1 =sqrt( std::pow(r1_hat,2) + std::pow(beta_new,2) );
00104                 c_old = c; // store for next iteration
00105                 s_old = s; // store for next iteration
00106                 c=r1_hat/r1; // new cosine
00107                 s=beta_new/r1; // new sine
00108                 
00109                 // Update solution
00110 //                p_oold = p_old;
00111                 const VectorType p_oold(p_old); // NOT SURE IF CREATING p_oold EVERY ITERATION IS EFFICIENT
00112                 p_old = p;
00113                 p.noalias()=(w-r2*p_old-r3*p_oold) /r1; // IS NOALIAS REQUIRED?
00114                 x += beta_one*c*eta*p;
00115                 residualNorm2 *= s*s;
00116                 
00117                 if ( residualNorm2 < threshold2){
00118                     break;
00119                 }
00120                 
00121                 eta=-s*eta; // update eta
00122                 iters++; // increment iteration number (for output purposes)
00123             }
00124             tol_error = std::sqrt(residualNorm2 / rhsNorm2); // return error. Note that this is the estimated error. The real error |Ax-b|/|b| may be slightly larger 
00125         }
00126         
00127     }
00128     
00129     template< typename _MatrixType, int _UpLo=Lower,
00130     typename _Preconditioner = IdentityPreconditioner>
00131 //    typename _Preconditioner = IdentityPreconditioner<typename _MatrixType::Scalar> > // preconditioner must be positive definite
00132     class MINRES;
00133     
00134     namespace internal {
00135         
00136         template< typename _MatrixType, int _UpLo, typename _Preconditioner>
00137         struct traits<MINRES<_MatrixType,_UpLo,_Preconditioner> >
00138         {
00139             typedef _MatrixType MatrixType;
00140             typedef _Preconditioner Preconditioner;
00141         };
00142         
00143     }
00144     
00194     template< typename _MatrixType, int _UpLo, typename _Preconditioner>
00195     class MINRES : public IterativeSolverBase<MINRES<_MatrixType,_UpLo,_Preconditioner> >
00196     {
00197         
00198         typedef IterativeSolverBase<MINRES> Base;
00199         using Base::mp_matrix;
00200         using Base::m_error;
00201         using Base::m_iterations;
00202         using Base::m_info;
00203         using Base::m_isInitialized;
00204     public:
00205         typedef _MatrixType MatrixType;
00206         typedef typename MatrixType::Scalar Scalar;
00207         typedef typename MatrixType::Index Index;
00208         typedef typename MatrixType::RealScalar RealScalar;
00209         typedef _Preconditioner Preconditioner;
00210         
00211         enum {UpLo = _UpLo};
00212         
00213     public:
00214         
00216         MINRES() : Base() {}
00217         
00228         MINRES(const MatrixType& A) : Base(A) {}
00229         
00231         ~MINRES(){}
00232                 
00238         template<typename Rhs,typename Guess>
00239         inline const internal::solve_retval_with_guess<MINRES, Rhs, Guess>
00240         solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
00241         {
00242             eigen_assert(m_isInitialized && "MINRES is not initialized.");
00243             eigen_assert(Base::rows()==b.rows()
00244                          && "MINRES::solve(): invalid number of rows of the right hand side matrix b");
00245             return internal::solve_retval_with_guess
00246             <MINRES, Rhs, Guess>(*this, b.derived(), x0);
00247         }
00248         
00250         template<typename Rhs,typename Dest>
00251         void _solveWithGuess(const Rhs& b, Dest& x) const
00252         {
00253             m_iterations = Base::maxIterations();
00254             m_error = Base::m_tolerance;
00255             
00256             for(int j=0; j<b.cols(); ++j)
00257             {
00258                 m_iterations = Base::maxIterations();
00259                 m_error = Base::m_tolerance;
00260                 
00261                 typename Dest::ColXpr xj(x,j);
00262                 internal::minres(mp_matrix->template selfadjointView<UpLo>(), b.col(j), xj,
00263                                  Base::m_preconditioner, m_iterations, m_error);
00264             }
00265             
00266             m_isInitialized = true;
00267             m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
00268         }
00269         
00271         template<typename Rhs,typename Dest>
00272         void _solve(const Rhs& b, Dest& x) const
00273         {
00274             x.setZero();
00275             _solveWithGuess(b,x);
00276         }
00277         
00278     protected:
00279         
00280     };
00281     
00282     namespace internal {
00283         
00284         template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs>
00285         struct solve_retval<MINRES<_MatrixType,_UpLo,_Preconditioner>, Rhs>
00286         : solve_retval_base<MINRES<_MatrixType,_UpLo,_Preconditioner>, Rhs>
00287         {
00288             typedef MINRES<_MatrixType,_UpLo,_Preconditioner> Dec;
00289             EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00290             
00291             template<typename Dest> void evalTo(Dest& dst) const
00292             {
00293                 dec()._solve(rhs(),dst);
00294             }
00295         };
00296         
00297     } // end namespace internal
00298     
00299 } // end namespace Eigen
00300 
00301 #endif // EIGEN_MINRES_H
00302 


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