12 #ifndef EIGEN_MINRES_H_
13 #define EIGEN_MINRES_H_
29 template<
typename MatrixType,
typename Rhs,
typename Dest,
typename Preconditioner>
31 void minres(
const MatrixType& mat,
const Rhs& rhs, Dest& x,
32 const Preconditioner& precond, Index& iters,
33 typename Dest::RealScalar& tol_error)
36 typedef typename Dest::RealScalar RealScalar;
37 typedef typename Dest::Scalar Scalar;
38 typedef Matrix<Scalar,Dynamic,1> VectorType;
41 const RealScalar rhsNorm2(rhs.squaredNorm());
51 const Index maxIters(iters);
52 const Index N(mat.cols());
53 const RealScalar threshold2(tol_error*tol_error*rhsNorm2);
57 VectorType v( VectorType::Zero(N) );
58 VectorType v_new(rhs-mat*x);
59 RealScalar residualNorm2(v_new.squaredNorm());
61 VectorType w_new(precond.solve(v_new));
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);
71 RealScalar c_old(1.0);
73 RealScalar s_old(0.0);
75 VectorType p_old(VectorType::Zero(N));
80 while ( iters < maxIters )
92 const RealScalar beta(beta_new);
98 v_new.noalias() = mat*w - beta*v_old;
99 const RealScalar alpha = v_new.dot(w);
101 w_new = precond.solve(v_new);
102 beta_new2 = v_new.dot(w_new);
103 eigen_assert(beta_new2 >= 0.0 &&
"PRECONDITIONER IS NOT POSITIVE DEFINITE");
104 beta_new = sqrt(beta_new2);
109 const RealScalar r2 =s*alpha+c*c_old*beta;
110 const RealScalar r3 =s_old*beta;
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) );
122 p.noalias()=(w-r2*p_old-r3*p_oold) /r1;
123 x += beta_one*c*eta*p;
127 residualNorm2 *= s*s;
129 if ( residualNorm2 < threshold2)
140 tol_error = std::sqrt(residualNorm2 / rhsNorm2);
145 template<
typename _MatrixType,
int _UpLo=Lower,
146 typename _Preconditioner = IdentityPreconditioner>
151 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
152 struct traits<
MINRES<_MatrixType,_UpLo,_Preconditioner> >
154 typedef _MatrixType MatrixType;
155 typedef _Preconditioner Preconditioner;
196 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
197 class MINRES :
public IterativeSolverBase<MINRES<_MatrixType,_UpLo,_Preconditioner> >
200 typedef IterativeSolverBase<MINRES> Base;
201 using Base::mp_matrix;
203 using Base::m_iterations;
205 using Base::m_isInitialized;
207 using Base::_solve_impl;
208 typedef _MatrixType MatrixType;
209 typedef typename MatrixType::Scalar Scalar;
210 typedef typename MatrixType::RealScalar RealScalar;
211 typedef _Preconditioner Preconditioner;
236 template<
typename Rhs,
typename Dest>
237 void _solve_with_guess_impl(
const Rhs& b, Dest& x)
const
239 typedef typename internal::conditional<UpLo==(Lower|Upper),
240 Ref<const MatrixType>&,
241 SparseSelfAdjointView<
const Ref<const MatrixType>, UpLo>
242 >::type MatrixWrapperType;
244 m_iterations = Base::maxIterations();
245 m_error = Base::m_tolerance;
247 for(
int j=0; j<b.cols(); ++j)
249 m_iterations = Base::maxIterations();
250 m_error = Base::m_tolerance;
252 typename Dest::ColXpr xj(x,j);
253 internal::minres(MatrixWrapperType(mp_matrix), b.col(j), xj,
254 Base::m_preconditioner, m_iterations, m_error);
257 m_isInitialized =
true;
258 m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
262 template<
typename Rhs,
typename Dest>
263 void _solve_impl(
const Rhs& b, MatrixBase<Dest> &x)
const
266 _solve_with_guess_impl(b,x.derived());
275 #endif // EIGEN_MINRES_H
Namespace containing all symbols from the Eigen library.
Definition: CXX11Meta.h:13
MINRES()
Definition: MINRES.h:218
A minimal residual solver for sparse symmetric problems.
Definition: MINRES.h:147
MINRES(const MatrixType &A)
Definition: MINRES.h:230
~MINRES()
Definition: MINRES.h:233