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;
198 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
199 class MINRES :
public IterativeSolverBase<MINRES<_MatrixType,_UpLo,_Preconditioner> >
202 typedef IterativeSolverBase<MINRES> Base;
205 using Base::m_iterations;
207 using Base::m_isInitialized;
209 using Base::_solve_impl;
210 typedef _MatrixType MatrixType;
211 typedef typename MatrixType::Scalar Scalar;
212 typedef typename MatrixType::RealScalar RealScalar;
213 typedef _Preconditioner Preconditioner;
232 template<
typename MatrixDerived>
233 explicit MINRES(
const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
239 template<
typename Rhs,
typename Dest>
240 void _solve_with_guess_impl(
const Rhs& b, Dest& x)
const 242 typedef typename Base::MatrixWrapper MatrixWrapper;
243 typedef typename Base::ActualMatrixType ActualMatrixType;
245 TransposeInput = (!MatrixWrapper::MatrixFree)
246 && (UpLo==(Lower|Upper))
247 && (!MatrixType::IsRowMajor)
248 && (!NumTraits<Scalar>::IsComplex)
250 typedef typename internal::conditional<TransposeInput,Transpose<const ActualMatrixType>, ActualMatrixType
const&>::type RowMajorWrapper;
251 EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(MatrixWrapper::MatrixFree,UpLo==(Lower|Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY);
252 typedef typename internal::conditional<UpLo==(Lower|Upper),
254 typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
255 >::type SelfAdjointWrapper;
257 m_iterations = Base::maxIterations();
258 m_error = Base::m_tolerance;
259 RowMajorWrapper row_mat(matrix());
260 for(
int j=0; j<b.cols(); ++j)
262 m_iterations = Base::maxIterations();
263 m_error = Base::m_tolerance;
265 typename Dest::ColXpr xj(x,j);
266 internal::minres(SelfAdjointWrapper(row_mat), b.col(j), xj,
267 Base::m_preconditioner, m_iterations, m_error);
270 m_isInitialized =
true;
271 m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
275 template<
typename Rhs,
typename Dest>
276 void _solve_impl(
const Rhs& b, MatrixBase<Dest> &x)
const 279 _solve_with_guess_impl(b,x.derived());
288 #endif // EIGEN_MINRES_H Namespace containing all symbols from the Eigen library.
Definition: AdolcForward:45
MINRES(const EigenBase< MatrixDerived > &A)
Definition: MINRES.h:233
MINRES()
Definition: MINRES.h:220
A minimal residual solver for sparse symmetric problems.
Definition: MINRES.h:147
~MINRES()
Definition: MINRES.h:236