35 template<
typename Vector,
typename RealScalar>
40 const RealScalar ns = s.norm();
41 const RealScalar nt = t.norm();
42 const Scalar ts = t.dot(s);
43 const RealScalar rho =
abs(ts / (nt * ns));
46 if (ts == Scalar(0)) {
53 return angle * (ns / nt) * (ts /
abs(ts));
55 return ts / (nt * nt);
58 template <
typename MatrixType,
typename Rhs,
typename Dest,
typename Preconditioner>
59 bool idrs(
const MatrixType& A,
const Rhs& b, Dest& x,
const Preconditioner& precond,
61 typename Dest::RealScalar& relres,
Index S,
bool smoothing,
typename Dest::RealScalar angle,
bool replacement)
63 typedef typename Dest::RealScalar RealScalar;
64 typedef typename Dest::Scalar Scalar;
66 typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> DenseMatrixType;
67 const Index N = b.size();
68 S = S < x.rows() ? S : x.rows();
69 const RealScalar tol = relres;
70 const Index maxit = iter;
72 Index replacements = 0;
75 FullPivLU<DenseMatrixType> lu_solver;
79 HouseholderQR<DenseMatrixType> qr(DenseMatrixType::Random(N, S));
80 P = (qr.householderQ() * DenseMatrixType::Identity(N, S));
83 const RealScalar normb = b.norm();
85 if (internal::isApprox(normb, RealScalar(0)))
104 const RealScalar mp = RealScalar(1e3) * NumTraits<Scalar>::epsilon();
109 const RealScalar tolb = tol * normb;
110 VectorType r = b - A * x;
120 RealScalar normr = r.norm();
126 relres = normr / normb;
130 DenseMatrixType G = DenseMatrixType::Zero(N, S);
131 DenseMatrixType U = DenseMatrixType::Zero(N, S);
132 DenseMatrixType M = DenseMatrixType::Identity(S, S);
133 VectorType t(N), v(N);
139 while (normr > tolb && iter < maxit)
142 VectorType f = (r.adjoint() * P).adjoint();
144 for (
Index k = 0; k < S; ++k)
148 lu_solver.compute(M.block(k , k , S -k, S - k ));
149 VectorType c = lu_solver.solve(f.segment(k , S - k ));
151 v = r - G.rightCols(S - k ) * c;
153 v = precond.solve(v);
156 U.col(k) = U.rightCols(S - k ) * c + om * v;
157 G.col(k) = A * U.col(k );
160 for (
Index i = 0; i < k-1 ; ++i)
163 Scalar alpha = P.col(i ).dot(G.col(k )) / M(i, i );
164 G.col(k ) = G.col(k ) - alpha * G.col(i );
165 U.col(k ) = U.col(k ) - alpha * U.col(i );
170 M.block(k , k , S - k , 1) = (G.col(k ).adjoint() * P.rightCols(S - k )).adjoint();
172 if (internal::isApprox(M(k,k), Scalar(0)))
178 Scalar beta = f(k ) / M(k , k );
179 r = r - beta * G.col(k );
180 x = x + beta * U.col(k );
183 if (replacement && normr > tolb / mp)
193 Scalar gamma = t.dot(r_s) / t.norm();
194 r_s = r_s - gamma * t;
195 x_s = x_s - gamma * (x_s - x);
199 if (normr < tolb || iter == maxit)
207 f.segment(k + 1, S - (k + 1) ) = f.segment(k + 1 , S - (k + 1)) - beta * M.block(k + 1 , k , S - (k + 1), 1);
211 if (normr < tolb || iter == maxit)
220 v = precond.solve(v);
226 om = internal::omega(t, r, angle);
228 if (om == RealScalar(0.0))
237 if (replacement && normr > tolb / mp)
243 if (trueres && normr < normb)
254 Scalar gamma = t.dot(r_s) /t.norm();
255 r_s = r_s - gamma * t;
256 x_s = x_s - gamma * (x_s - x);
274 template <
typename _MatrixType,
typename _Preconditioner = DiagonalPreconditioner<
typename _MatrixType::Scalar> >
280 template <
typename _MatrixType,
typename _Preconditioner>
281 struct traits<
Eigen::
IDRS<_MatrixType, _Preconditioner> >
283 typedef _MatrixType MatrixType;
284 typedef _Preconditioner Preconditioner;
330 template <
typename _MatrixType,
typename _Preconditioner>
335 typedef _MatrixType MatrixType;
336 typedef typename MatrixType::Scalar Scalar;
337 typedef typename MatrixType::RealScalar RealScalar;
338 typedef _Preconditioner Preconditioner;
344 using Base::m_isInitialized;
345 using Base::m_iterations;
354 IDRS(): m_S(4), m_smoothing(false), m_angle(RealScalar(0.7)), m_residual(false) {}
366 template <
typename MatrixDerived>
368 m_angle(RealScalar(0.7)), m_residual(false) {}
376 template <
typename Rhs,
typename Dest>
379 m_iterations = Base::maxIterations();
380 m_error = Base::m_tolerance;
382 bool ret = internal::idrs(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error, m_S,m_smoothing,m_angle,m_residual);
406 m_smoothing=smoothing;
void setSmoothing(bool smoothing)
Definition: IDRS.h:404
internal::traits< Derived >::Scalar Scalar
IDRS()
Definition: IDRS.h:354
Namespace containing all symbols from the Eigen library.
The Induced Dimension Reduction method (IDR(s)) is a short-recurrences Krylov method for sparse squar...
Definition: IDRS.h:275
Matrix< Type, Size, 1 > Vector
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
void _solve_vector_with_guess_impl(const Rhs &b, Dest &x) const
Definition: IDRS.h:377
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
void setResidualUpdate(bool update)
Definition: IDRS.h:427
void setS(Index S)
Definition: IDRS.h:388
IDRS(const EigenBase< MatrixDerived > &A)
Definition: IDRS.h:367
void setAngle(RealScalar angle)
Definition: IDRS.h:419