13 #ifndef EIGEN_MINRES_H_ 14 #define EIGEN_MINRES_H_ 30 template<
typename MatrixType,
typename Rhs,
typename Dest,
typename Preconditioner>
32 void minres(
const MatrixType& mat,
const Rhs& rhs, Dest& x,
33 const Preconditioner& precond,
Index& iters,
34 typename Dest::RealScalar& tol_error)
37 typedef typename Dest::RealScalar RealScalar;
38 typedef typename Dest::Scalar Scalar;
42 const RealScalar rhsNorm2(rhs.squaredNorm());
52 const Index maxIters(iters);
53 const Index N(mat.cols());
54 const RealScalar threshold2(tol_error*tol_error*rhsNorm2);
58 VectorType v( VectorType::Zero(N) );
59 VectorType v_new(rhs-mat*x);
60 RealScalar residualNorm2(v_new.squaredNorm());
62 VectorType w_new(precond.solve(v_new));
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);
70 RealScalar c_old(1.0);
72 RealScalar s_old(0.0);
74 VectorType p_old(VectorType::Zero(N));
79 while ( iters < maxIters )
91 const RealScalar beta(beta_new);
97 v_new.noalias() = mat*w - beta*v_old;
98 const RealScalar alpha = v_new.dot(w);
100 w_new = precond.solve(v_new);
101 beta_new2 = v_new.dot(w_new);
102 eigen_assert(beta_new2 >= 0.0 &&
"PRECONDITIONER IS NOT POSITIVE DEFINITE");
103 beta_new =
sqrt(beta_new2);
106 const RealScalar r2 =s*alpha+c*c_old*beta;
107 const RealScalar r3 =s_old*beta;
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) );
118 p.noalias()=(w-r2*p_old-r3*p_oold) /r1;
119 x += beta_one*c*eta*p;
123 residualNorm2 *= s*s;
125 if ( residualNorm2 < threshold2)
136 tol_error = std::sqrt(residualNorm2 / rhsNorm2);
141 template<
typename _MatrixType,
int _UpLo=
Lower,
142 typename _Preconditioner = IdentityPreconditioner>
147 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
148 struct traits<
MINRES<_MatrixType,_UpLo,_Preconditioner> >
150 typedef _MatrixType MatrixType;
151 typedef _Preconditioner Preconditioner;
194 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
201 using Base::m_iterations;
203 using Base::m_isInitialized;
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;
228 template<
typename MatrixDerived>
235 template<
typename Rhs,
typename Dest>
236 void _solve_vector_with_guess_impl(
const Rhs& b, Dest& x)
const 239 typedef typename Base::ActualMatrixType ActualMatrixType;
241 TransposeInput = (!MatrixWrapper::MatrixFree)
243 && (!MatrixType::IsRowMajor)
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),
250 typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
251 >::type SelfAdjointWrapper;
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);
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