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,
int& 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 int maxIters(iters);
42 const int N(mat.cols());
43 const RealScalar rhsNorm2(rhs.squaredNorm());
44 const RealScalar threshold2(tol_error*tol_error*rhsNorm2);
48 VectorType v( VectorType::Zero(N) );
49 VectorType v_new(rhs-mat*x);
50 RealScalar residualNorm2(v_new.squaredNorm());
52 VectorType w_new(precond.solve(v_new));
54 RealScalar beta_new2(v_new.dot(w_new));
55 eigen_assert(beta_new2 >= 0 &&
"PRECONDITIONER IS NOT POSITIVE DEFINITE");
56 RealScalar beta_new(sqrt(beta_new2));
57 const RealScalar beta_one(beta_new);
62 RealScalar c_old(1.0);
64 RealScalar s_old(0.0);
66 VectorType p_old(VectorType::Zero(N));
71 while ( iters < maxIters ){
83 const RealScalar beta(beta_new);
85 const VectorType v_old(v);
88 const VectorType w(w_new);
89 v_new.noalias() = mat*w - beta*v_old;
90 const RealScalar alpha = v_new.dot(w);
92 w_new = precond.solve(v_new);
93 beta_new2 = v_new.dot(w_new);
94 eigen_assert(beta_new2 >= 0 &&
"PRECONDITIONER IS NOT POSITIVE DEFINITE");
95 beta_new = sqrt(beta_new2);
100 const RealScalar r2 =s*alpha+c*c_old*beta;
101 const RealScalar r3 =s_old*beta;
102 const RealScalar r1_hat=c*alpha-c_old*s*beta;
103 const RealScalar r1 =sqrt( std::pow(r1_hat,2) + std::pow(beta_new,2) );
111 const VectorType p_oold(p_old);
113 p.noalias()=(w-r2*p_old-r3*p_oold) /r1;
114 x += beta_one*c*eta*p;
115 residualNorm2 *= s*s;
117 if ( residualNorm2 < threshold2){
124 tol_error = std::sqrt(residualNorm2 / rhsNorm2);
129 template<
typename _MatrixType,
int _UpLo=Lower,
130 typename _Preconditioner = IdentityPreconditioner>
136 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
137 struct traits<
MINRES<_MatrixType,_UpLo,_Preconditioner> >
139 typedef _MatrixType MatrixType;
140 typedef _Preconditioner Preconditioner;
194 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
195 class MINRES :
public IterativeSolverBase<MINRES<_MatrixType,_UpLo,_Preconditioner> >
198 typedef IterativeSolverBase<MINRES> Base;
199 using Base::mp_matrix;
201 using Base::m_iterations;
203 using Base::m_isInitialized;
205 typedef _MatrixType MatrixType;
206 typedef typename MatrixType::Scalar Scalar;
207 typedef typename MatrixType::Index Index;
208 typedef typename MatrixType::RealScalar RealScalar;
209 typedef _Preconditioner Preconditioner;
238 template<
typename Rhs,
typename Guess>
239 inline const internal::solve_retval_with_guess<MINRES, Rhs, Guess>
242 eigen_assert(m_isInitialized &&
"MINRES is not initialized.");
243 eigen_assert(Base::rows()==b.rows()
244 &&
"MINRES::solve(): invalid number of rows of the right hand side matrix b");
245 return internal::solve_retval_with_guess
246 <
MINRES, Rhs, Guess>(*
this, b.derived(), x0);
250 template<
typename Rhs,
typename Dest>
251 void _solveWithGuess(
const Rhs& b, Dest& x)
const 253 m_iterations = Base::maxIterations();
254 m_error = Base::m_tolerance;
256 for(
int j=0; j<b.cols(); ++j)
258 m_iterations = Base::maxIterations();
259 m_error = Base::m_tolerance;
261 typename Dest::ColXpr xj(x,j);
262 internal::minres(mp_matrix->template selfadjointView<UpLo>(), b.col(j), xj,
263 Base::m_preconditioner, m_iterations, m_error);
266 m_isInitialized =
true;
267 m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
271 template<
typename Rhs,
typename Dest>
272 void _solve(
const Rhs& b, Dest& x)
const 275 _solveWithGuess(b,x);
284 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner,
typename Rhs>
285 struct solve_retval<MINRES<_MatrixType,_UpLo,_Preconditioner>, Rhs>
286 : solve_retval_base<MINRES<_MatrixType,_UpLo,_Preconditioner>, Rhs>
289 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
291 template<
typename Dest>
void evalTo(Dest& dst)
const 293 dec()._solve(rhs(),dst);
301 #endif // EIGEN_MINRES_H const internal::solve_retval_with_guess< MINRES, Rhs, Guess > solveWithGuess(const MatrixBase< Rhs > &b, const Guess &x0) const
Definition: MINRES.h:240
~MINRES()
Definition: MINRES.h:231
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: AutoDiffJacobian.h:13
A minimal residual solver for sparse symmetric problems.
Definition: MINRES.h:132
MINRES()
Definition: MINRES.h:216
MINRES(const MatrixType &A)
Definition: MINRES.h:228