Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
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-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
6 // Copyright (C) 2018 David Hyde <dabh@stanford.edu>
7 //
8 // This Source Code Form is subject to the terms of the Mozilla
9 // Public License v. 2.0. If a copy of the MPL was not distributed
10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
11 
12 
13 #ifndef EIGEN_MINRES_H_
14 #define EIGEN_MINRES_H_
15 
16 
17 namespace Eigen {
18 
19  namespace internal {
20 
30  template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
31  EIGEN_DONT_INLINE
32  void minres(const MatrixType& mat, const Rhs& rhs, Dest& x,
33  const Preconditioner& precond, Index& iters,
34  typename Dest::RealScalar& tol_error)
35  {
36  using std::sqrt;
37  typedef typename Dest::RealScalar RealScalar;
38  typedef typename Dest::Scalar Scalar;
39  typedef Matrix<Scalar,Dynamic,1> VectorType;
40 
41  // Check for zero rhs
42  const RealScalar rhsNorm2(rhs.squaredNorm());
43  if(rhsNorm2 == 0)
44  {
45  x.setZero();
46  iters = 0;
47  tol_error = 0;
48  return;
49  }
50 
51  // initialize
52  const Index maxIters(iters); // initialize maxIters to iters
53  const Index N(mat.cols()); // the size of the matrix
54  const RealScalar threshold2(tol_error*tol_error*rhsNorm2); // convergence threshold (compared to residualNorm2)
55 
56  // Initialize preconditioned Lanczos
57  VectorType v_old(N); // will be initialized inside loop
58  VectorType v( VectorType::Zero(N) ); //initialize v
59  VectorType v_new(rhs-mat*x); //initialize v_new
60  RealScalar residualNorm2(v_new.squaredNorm());
61  VectorType w(N); // will be initialized inside loop
62  VectorType w_new(precond.solve(v_new)); // initialize w_new
63 // RealScalar beta; // will be initialized inside loop
64  RealScalar beta_new2(v_new.dot(w_new));
65  eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
66  RealScalar beta_new(sqrt(beta_new2));
67  const RealScalar beta_one(beta_new);
68  // Initialize other variables
69  RealScalar c(1.0); // the cosine of the Givens rotation
70  RealScalar c_old(1.0);
71  RealScalar s(0.0); // the sine of the Givens rotation
72  RealScalar s_old(0.0); // the sine of the Givens rotation
73  VectorType p_oold(N); // will be initialized in loop
74  VectorType p_old(VectorType::Zero(N)); // initialize p_old=0
75  VectorType p(p_old); // initialize p=0
76  RealScalar eta(1.0);
77 
78  iters = 0; // reset iters
79  while ( iters < maxIters )
80  {
81  // Preconditioned Lanczos
82  /* Note that there are 4 variants on the Lanczos algorithm. These are
83  * described in Paige, C. C. (1972). Computational variants of
84  * the Lanczos method for the eigenproblem. IMA Journal of Applied
85  * Mathematics, 10(3), 373-381. The current implementation corresponds
86  * to the case A(2,7) in the paper. It also corresponds to
87  * algorithm 6.14 in Y. Saad, Iterative Methods for Sparse Linear
88  * Systems, 2003 p.173. For the preconditioned version see
89  * A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM (1987).
90  */
91  const RealScalar beta(beta_new);
92  v_old = v; // update: at first time step, this makes v_old = 0 so value of beta doesn't matter
93  v_new /= beta_new; // overwrite v_new for next iteration
94  w_new /= beta_new; // overwrite w_new for next iteration
95  v = v_new; // update
96  w = w_new; // update
97  v_new.noalias() = mat*w - beta*v_old; // compute v_new
98  const RealScalar alpha = v_new.dot(w);
99  v_new -= alpha*v; // overwrite v_new
100  w_new = precond.solve(v_new); // overwrite w_new
101  beta_new2 = v_new.dot(w_new); // compute beta_new
102  eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
103  beta_new = sqrt(beta_new2); // compute beta_new
104 
105  // Givens rotation
106  const RealScalar r2 =s*alpha+c*c_old*beta; // s, s_old, c and c_old are still from previous iteration
107  const RealScalar r3 =s_old*beta; // s, s_old, c and c_old are still from previous iteration
108  const RealScalar r1_hat=c*alpha-c_old*s*beta;
109  const RealScalar r1 =sqrt( std::pow(r1_hat,2) + std::pow(beta_new,2) );
110  c_old = c; // store for next iteration
111  s_old = s; // store for next iteration
112  c=r1_hat/r1; // new cosine
113  s=beta_new/r1; // new sine
114 
115  // Update solution
116  p_oold = p_old;
117  p_old = p;
118  p.noalias()=(w-r2*p_old-r3*p_oold) /r1; // IS NOALIAS REQUIRED?
119  x += beta_one*c*eta*p;
120 
121  /* Update the squared residual. Note that this is the estimated residual.
122  The real residual |Ax-b|^2 may be slightly larger */
123  residualNorm2 *= s*s;
124 
125  if ( residualNorm2 < threshold2)
126  {
127  break;
128  }
129 
130  eta=-s*eta; // update eta
131  iters++; // increment iteration number (for output purposes)
132  }
133 
134  /* Compute error. Note that this is the estimated error. The real
135  error |Ax-b|/|b| may be slightly larger */
136  tol_error = std::sqrt(residualNorm2 / rhsNorm2);
137  }
138 
139  }
140 
141  template< typename _MatrixType, int _UpLo=Lower,
142  typename _Preconditioner = IdentityPreconditioner>
143  class MINRES;
144 
145  namespace internal {
146 
147  template< typename _MatrixType, int _UpLo, typename _Preconditioner>
148  struct traits<MINRES<_MatrixType,_UpLo,_Preconditioner> >
149  {
150  typedef _MatrixType MatrixType;
151  typedef _Preconditioner Preconditioner;
152  };
153 
154  }
155 
194  template< typename _MatrixType, int _UpLo, typename _Preconditioner>
195  class MINRES : public IterativeSolverBase<MINRES<_MatrixType,_UpLo,_Preconditioner> >
196  {
197 
199  using Base::matrix;
200  using Base::m_error;
201  using Base::m_iterations;
202  using Base::m_info;
203  using Base::m_isInitialized;
204  public:
205  using Base::_solve_impl;
206  typedef _MatrixType MatrixType;
207  typedef typename MatrixType::Scalar Scalar;
208  typedef typename MatrixType::RealScalar RealScalar;
209  typedef _Preconditioner Preconditioner;
210 
211  enum {UpLo = _UpLo};
212 
213  public:
214 
216  MINRES() : Base() {}
217 
228  template<typename MatrixDerived>
229  explicit MINRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
230 
233 
235  template<typename Rhs,typename Dest>
236  void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
237  {
238  typedef typename Base::MatrixWrapper MatrixWrapper;
239  typedef typename Base::ActualMatrixType ActualMatrixType;
240  enum {
241  TransposeInput = (!MatrixWrapper::MatrixFree)
242  && (UpLo==(Lower|Upper))
243  && (!MatrixType::IsRowMajor)
245  };
246  typedef typename internal::conditional<TransposeInput,Transpose<const ActualMatrixType>, ActualMatrixType const&>::type RowMajorWrapper;
247  EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(MatrixWrapper::MatrixFree,UpLo==(Lower|Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY);
248  typedef typename internal::conditional<UpLo==(Lower|Upper),
249  RowMajorWrapper,
250  typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
251  >::type SelfAdjointWrapper;
252 
253  m_iterations = Base::maxIterations();
254  m_error = Base::m_tolerance;
255  RowMajorWrapper row_mat(matrix());
256  internal::minres(SelfAdjointWrapper(row_mat), b, x,
257  Base::m_preconditioner, m_iterations, m_error);
258  m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
259  }
260 
261  protected:
262 
263  };
264 
265 } // end namespace Eigen
266 
267 #endif // EIGEN_MINRES_H
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Namespace containing all symbols from the Eigen library.
MINRES(const EigenBase< MatrixDerived > &A)
Definition: MINRES.h:229
MINRES()
Definition: MINRES.h:216
A minimal residual solver for sparse symmetric problems.
Definition: MINRES.h:143
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
~MINRES()
Definition: MINRES.h:232