Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
 
Loading...
Searching...
No Matches
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
17namespace 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
A minimal residual solver for sparse symmetric problems.
Definition: MINRES.h:196
~MINRES()
Definition: MINRES.h:232
MINRES(const EigenBase< MatrixDerived > &A)
Definition: MINRES.h:229
MINRES()
Definition: MINRES.h:216
Namespace containing all symbols from the Eigen library.
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index