All Classes Namespaces Functions Variables Typedefs Enumerator Groups Pages
MINRES.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012 Giacomo Po <gpo@ucla.edu>
5 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 
12 #ifndef EIGEN_MINRES_H_
13 #define EIGEN_MINRES_H_
14 
15 
16 namespace Eigen {
17 
18  namespace internal {
19 
29  template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
30  EIGEN_DONT_INLINE
31  void minres(const MatrixType& mat, const Rhs& rhs, Dest& x,
32  const Preconditioner& precond, int& iters,
33  typename Dest::RealScalar& tol_error)
34  {
35  using std::sqrt;
36  typedef typename Dest::RealScalar RealScalar;
37  typedef typename Dest::Scalar Scalar;
38  typedef Matrix<Scalar,Dynamic,1> VectorType;
39 
40  // Check for zero rhs
41  const RealScalar rhsNorm2(rhs.squaredNorm());
42  if(rhsNorm2 == 0)
43  {
44  x.setZero();
45  iters = 0;
46  tol_error = 0;
47  return;
48  }
49 
50  // initialize
51  const int maxIters(iters); // initialize maxIters to iters
52  const int N(mat.cols()); // the size of the matrix
53  const RealScalar threshold2(tol_error*tol_error*rhsNorm2); // convergence threshold (compared to residualNorm2)
54 
55  // Initialize preconditioned Lanczos
56  VectorType v_old(N); // will be initialized inside loop
57  VectorType v( VectorType::Zero(N) ); //initialize v
58  VectorType v_new(rhs-mat*x); //initialize v_new
59  RealScalar residualNorm2(v_new.squaredNorm());
60  VectorType w(N); // will be initialized inside loop
61  VectorType w_new(precond.solve(v_new)); // initialize w_new
62 // RealScalar beta; // will be initialized inside loop
63  RealScalar beta_new2(v_new.dot(w_new));
64  eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
65  RealScalar beta_new(sqrt(beta_new2));
66  const RealScalar beta_one(beta_new);
67  v_new /= beta_new;
68  w_new /= beta_new;
69  // Initialize other variables
70  RealScalar c(1.0); // the cosine of the Givens rotation
71  RealScalar c_old(1.0);
72  RealScalar s(0.0); // the sine of the Givens rotation
73  RealScalar s_old(0.0); // the sine of the Givens rotation
74  VectorType p_oold(N); // will be initialized in loop
75  VectorType p_old(VectorType::Zero(N)); // initialize p_old=0
76  VectorType p(p_old); // initialize p=0
77  RealScalar eta(1.0);
78 
79  iters = 0; // reset iters
80  while ( iters < maxIters )
81  {
82  // Preconditioned Lanczos
83  /* Note that there are 4 variants on the Lanczos algorithm. These are
84  * described in Paige, C. C. (1972). Computational variants of
85  * the Lanczos method for the eigenproblem. IMA Journal of Applied
86  * Mathematics, 10(3), 373–381. The current implementation corresponds
87  * to the case A(2,7) in the paper. It also corresponds to
88  * algorithm 6.14 in Y. Saad, Iterative Methods for Sparse Linear
89  * Systems, 2003 p.173. For the preconditioned version see
90  * A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM (1987).
91  */
92  const RealScalar beta(beta_new);
93  v_old = v; // update: at first time step, this makes v_old = 0 so value of beta doesn't matter
94 // const VectorType v_old(v); // NOT SURE IF CREATING v_old EVERY ITERATION IS EFFICIENT
95  v = v_new; // update
96  w = w_new; // update
97 // const VectorType w(w_new); // NOT SURE IF CREATING w EVERY ITERATION IS EFFICIENT
98  v_new.noalias() = mat*w - beta*v_old; // compute v_new
99  const RealScalar alpha = v_new.dot(w);
100  v_new -= alpha*v; // overwrite v_new
101  w_new = precond.solve(v_new); // overwrite w_new
102  beta_new2 = v_new.dot(w_new); // compute beta_new
103  eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
104  beta_new = sqrt(beta_new2); // compute beta_new
105  v_new /= beta_new; // overwrite v_new for next iteration
106  w_new /= beta_new; // overwrite w_new for next iteration
107 
108  // Givens rotation
109  const RealScalar r2 =s*alpha+c*c_old*beta; // s, s_old, c and c_old are still from previous iteration
110  const RealScalar r3 =s_old*beta; // s, s_old, c and c_old are still from previous iteration
111  const RealScalar r1_hat=c*alpha-c_old*s*beta;
112  const RealScalar r1 =sqrt( std::pow(r1_hat,2) + std::pow(beta_new,2) );
113  c_old = c; // store for next iteration
114  s_old = s; // store for next iteration
115  c=r1_hat/r1; // new cosine
116  s=beta_new/r1; // new sine
117 
118  // Update solution
119  p_oold = p_old;
120 // const VectorType p_oold(p_old); // NOT SURE IF CREATING p_oold EVERY ITERATION IS EFFICIENT
121  p_old = p;
122  p.noalias()=(w-r2*p_old-r3*p_oold) /r1; // IS NOALIAS REQUIRED?
123  x += beta_one*c*eta*p;
124 
125  /* Update the squared residual. Note that this is the estimated residual.
126  The real residual |Ax-b|^2 may be slightly larger */
127  residualNorm2 *= s*s;
128 
129  if ( residualNorm2 < threshold2)
130  {
131  break;
132  }
133 
134  eta=-s*eta; // update eta
135  iters++; // increment iteration number (for output purposes)
136  }
137 
138  /* Compute error. Note that this is the estimated error. The real
139  error |Ax-b|/|b| may be slightly larger */
140  tol_error = std::sqrt(residualNorm2 / rhsNorm2);
141  }
142 
143  }
144 
145  template< typename _MatrixType, int _UpLo=Lower,
146  typename _Preconditioner = IdentityPreconditioner>
147 // typename _Preconditioner = IdentityPreconditioner<typename _MatrixType::Scalar> > // preconditioner must be positive definite
148  class MINRES;
149 
150  namespace internal {
151 
152  template< typename _MatrixType, int _UpLo, typename _Preconditioner>
153  struct traits<MINRES<_MatrixType,_UpLo,_Preconditioner> >
154  {
155  typedef _MatrixType MatrixType;
156  typedef _Preconditioner Preconditioner;
157  };
158 
159  }
160 
197  template< typename _MatrixType, int _UpLo, typename _Preconditioner>
198  class MINRES : public IterativeSolverBase<MINRES<_MatrixType,_UpLo,_Preconditioner> >
199  {
200 
201  typedef IterativeSolverBase<MINRES> Base;
202  using Base::mp_matrix;
203  using Base::m_error;
204  using Base::m_iterations;
205  using Base::m_info;
206  using Base::m_isInitialized;
207  public:
208  typedef _MatrixType MatrixType;
209  typedef typename MatrixType::Scalar Scalar;
210  typedef typename MatrixType::Index Index;
211  typedef typename MatrixType::RealScalar RealScalar;
212  typedef _Preconditioner Preconditioner;
213 
214  enum {UpLo = _UpLo};
215 
216  public:
217 
219  MINRES() : Base() {}
220 
231  MINRES(const MatrixType& A) : Base(A) {}
232 
235 
241  template<typename Rhs,typename Guess>
242  inline const internal::solve_retval_with_guess<MINRES, Rhs, Guess>
243  solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
244  {
245  eigen_assert(m_isInitialized && "MINRES is not initialized.");
246  eigen_assert(Base::rows()==b.rows()
247  && "MINRES::solve(): invalid number of rows of the right hand side matrix b");
248  return internal::solve_retval_with_guess
249  <MINRES, Rhs, Guess>(*this, b.derived(), x0);
250  }
251 
253  template<typename Rhs,typename Dest>
254  void _solveWithGuess(const Rhs& b, Dest& x) const
255  {
256  typedef typename internal::conditional<UpLo==(Lower|Upper),
257  const MatrixType&,
258  SparseSelfAdjointView<const MatrixType, UpLo>
259  >::type MatrixWrapperType;
260 
261  m_iterations = Base::maxIterations();
262  m_error = Base::m_tolerance;
263 
264  for(int j=0; j<b.cols(); ++j)
265  {
266  m_iterations = Base::maxIterations();
267  m_error = Base::m_tolerance;
268 
269  typename Dest::ColXpr xj(x,j);
270  internal::minres(MatrixWrapperType(*mp_matrix), b.col(j), xj,
271  Base::m_preconditioner, m_iterations, m_error);
272  }
273 
274  m_isInitialized = true;
275  m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
276  }
277 
279  template<typename Rhs,typename Dest>
280  void _solve(const Rhs& b, Dest& x) const
281  {
282  x.setZero();
283  _solveWithGuess(b,x);
284  }
285 
286  protected:
287 
288  };
289 
290  namespace internal {
291 
292  template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs>
293  struct solve_retval<MINRES<_MatrixType,_UpLo,_Preconditioner>, Rhs>
294  : solve_retval_base<MINRES<_MatrixType,_UpLo,_Preconditioner>, Rhs>
295  {
296  typedef MINRES<_MatrixType,_UpLo,_Preconditioner> Dec;
297  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
298 
299  template<typename Dest> void evalTo(Dest& dst) const
300  {
301  dec()._solve(rhs(),dst);
302  }
303  };
304 
305  } // end namespace internal
306 
307 } // end namespace Eigen
308 
309 #endif // EIGEN_MINRES_H
310 
const internal::solve_retval_with_guess< MINRES, Rhs, Guess > solveWithGuess(const MatrixBase< Rhs > &b, const Guess &x0) const
Definition: MINRES.h:243
~MINRES()
Definition: MINRES.h:234
A minimal residual solver for sparse symmetric problems.
Definition: MINRES.h:148
MINRES()
Definition: MINRES.h:219
MINRES(const MatrixType &A)
Definition: MINRES.h:231