Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
 
Loading...
Searching...
No Matches
IDRS.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2020 Chris Schoutrop <c.e.m.schoutrop@tue.nl>
5// Copyright (C) 2020 Jens Wehner <j.wehner@esciencecenter.nl>
6// Copyright (C) 2020 Jan van Dijk <j.v.dijk@tue.nl>
7//
8// This Source Code Form is subject to the terms of the Mozilla
9// Public License v. 2.0. If a copy of the MPL was not distributed
10// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
11
12
13#ifndef EIGEN_IDRS_H
14#define EIGEN_IDRS_H
15
16namespace Eigen
17{
18
19 namespace internal
20 {
35 template<typename Vector, typename RealScalar>
36 typename Vector::Scalar omega(const Vector& t, const Vector& s, RealScalar angle)
37 {
38 using numext::abs;
39 typedef typename Vector::Scalar Scalar;
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));
44
45 if (rho < angle) {
46 if (ts == Scalar(0)) {
47 return Scalar(0);
48 }
49 // Original relation for om is given by
50 // om = om * angle / rho;
51 // To alleviate potential (near) division by zero this can be rewritten as
52 // om = angle * (ns / nt) * (ts / abs(ts)) = angle * (ns / nt) * sgn(ts)
53 return angle * (ns / nt) * (ts / abs(ts));
54 }
55 return ts / (nt * nt);
56 }
57
58 template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
59 bool idrs(const MatrixType& A, const Rhs& b, Dest& x, const Preconditioner& precond,
60 Index& iter,
61 typename Dest::RealScalar& relres, Index S, bool smoothing, typename Dest::RealScalar angle, bool replacement)
62 {
63 typedef typename Dest::RealScalar RealScalar;
64 typedef typename Dest::Scalar Scalar;
65 typedef Matrix<Scalar, Dynamic, 1> VectorType;
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;
71
72 Index replacements = 0;
73 bool trueres = false;
74
75 FullPivLU<DenseMatrixType> lu_solver;
76
77 DenseMatrixType P;
78 {
79 HouseholderQR<DenseMatrixType> qr(DenseMatrixType::Random(N, S));
80 P = (qr.householderQ() * DenseMatrixType::Identity(N, S));
81 }
82
83 const RealScalar normb = b.norm();
84
85 if (internal::isApprox(normb, RealScalar(0)))
86 {
87 //Solution is the zero vector
88 x.setZero();
89 iter = 0;
90 relres = 0;
91 return true;
92 }
93 // from http://homepage.tudelft.nl/1w5b5/IDRS/manual.pdf
94 // A peak in the residual is considered dangerously high if‖ri‖/‖b‖> C(tol/epsilon).
95 // With epsilon the
96 // relative machine precision. The factor tol/epsilon corresponds to the size of a
97 // finite precision number that is so large that the absolute round-off error in
98 // this number, when propagated through the process, makes it impossible to
99 // achieve the required accuracy.The factor C accounts for the accumulation of
100 // round-off errors. This parameter has beenset to 10−3.
101 // mp is epsilon/C
102 // 10^3 * eps is very conservative, so normally no residual replacements will take place.
103 // It only happens if things go very wrong. Too many restarts may ruin the convergence.
104 const RealScalar mp = RealScalar(1e3) * NumTraits<Scalar>::epsilon();
105
106
107
108 //Compute initial residual
109 const RealScalar tolb = tol * normb; //Relative tolerance
110 VectorType r = b - A * x;
111
112 VectorType x_s, r_s;
113
114 if (smoothing)
115 {
116 x_s = x;
117 r_s = r;
118 }
119
120 RealScalar normr = r.norm();
121
122 if (normr <= tolb)
123 {
124 //Initial guess is a good enough solution
125 iter = 0;
126 relres = normr / normb;
127 return true;
128 }
129
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);
134 Scalar om = 1.;
135
136 //Main iteration loop, guild G-spaces:
137 iter = 0;
138
139 while (normr > tolb && iter < maxit)
140 {
141 //New right hand size for small system:
142 VectorType f = (r.adjoint() * P).adjoint();
143
144 for (Index k = 0; k < S; ++k)
145 {
146 //Solve small system and make v orthogonal to P:
147 //c = M(k:s,k:s)\f(k:s);
148 lu_solver.compute(M.block(k , k , S -k, S - k ));
149 VectorType c = lu_solver.solve(f.segment(k , S - k ));
150 //v = r - G(:,k:s)*c;
151 v = r - G.rightCols(S - k ) * c;
152 //Preconditioning
153 v = precond.solve(v);
154
155 //Compute new U(:,k) and G(:,k), G(:,k) is in space G_j
156 U.col(k) = U.rightCols(S - k ) * c + om * v;
157 G.col(k) = A * U.col(k );
158
159 //Bi-Orthogonalise the new basis vectors:
160 for (Index i = 0; i < k-1 ; ++i)
161 {
162 //alpha = ( P(:,i)'*G(:,k) )/M(i,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 );
166 }
167
168 //New column of M = P'*G (first k-1 entries are zero)
169 //M(k:s,k) = (G(:,k)'*P(:,k:s))';
170 M.block(k , k , S - k , 1) = (G.col(k ).adjoint() * P.rightCols(S - k )).adjoint();
171
172 if (internal::isApprox(M(k,k), Scalar(0)))
173 {
174 return false;
175 }
176
177 //Make r orthogonal to q_i, i = 0..k-1
178 Scalar beta = f(k ) / M(k , k );
179 r = r - beta * G.col(k );
180 x = x + beta * U.col(k );
181 normr = r.norm();
182
183 if (replacement && normr > tolb / mp)
184 {
185 trueres = true;
186 }
187
188 //Smoothing:
189 if (smoothing)
190 {
191 t = r_s - r;
192 //gamma is a Scalar, but the conversion is not allowed
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);
196 normr = r_s.norm();
197 }
198
199 if (normr < tolb || iter == maxit)
200 {
201 break;
202 }
203
204 //New f = P'*r (first k components are zero)
205 if (k < S-1)
206 {
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);
208 }
209 }//end for
210
211 if (normr < tolb || iter == maxit)
212 {
213 break;
214 }
215
216 //Now we have sufficient vectors in G_j to compute residual in G_j+1
217 //Note: r is already perpendicular to P so v = r
218 //Preconditioning
219 v = r;
220 v = precond.solve(v);
221
222 //Matrix-vector multiplication:
223 t = A * v;
224
225 //Computation of a new omega
226 om = internal::omega(t, r, angle);
227
228 if (om == RealScalar(0.0))
229 {
230 return false;
231 }
232
233 r = r - om * t;
234 x = x + om * v;
235 normr = r.norm();
236
237 if (replacement && normr > tolb / mp)
238 {
239 trueres = true;
240 }
241
242 //Residual replacement?
243 if (trueres && normr < normb)
244 {
245 r = b - A * x;
246 trueres = false;
247 replacements++;
248 }
249
250 //Smoothing:
251 if (smoothing)
252 {
253 t = r_s - r;
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);
257 normr = r_s.norm();
258 }
259
260 iter++;
261
262 }//end while
263
264 if (smoothing)
265 {
266 x = x_s;
267 }
268 relres=normr/normb;
269 return true;
270 }
271
272 } // namespace internal
273
274 template <typename _MatrixType, typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
275 class IDRS;
276
277 namespace internal
278 {
279
280 template <typename _MatrixType, typename _Preconditioner>
281 struct traits<Eigen::IDRS<_MatrixType, _Preconditioner> >
282 {
283 typedef _MatrixType MatrixType;
284 typedef _Preconditioner Preconditioner;
285 };
286
287 } // namespace internal
288
289
330 template <typename _MatrixType, typename _Preconditioner>
331 class IDRS : public IterativeSolverBase<IDRS<_MatrixType, _Preconditioner> >
332 {
333
334 public:
335 typedef _MatrixType MatrixType;
336 typedef typename MatrixType::Scalar Scalar;
337 typedef typename MatrixType::RealScalar RealScalar;
338 typedef _Preconditioner Preconditioner;
339
340 private:
342 using Base::m_error;
343 using Base::m_info;
344 using Base::m_isInitialized;
345 using Base::m_iterations;
346 using Base::matrix;
347 Index m_S;
348 bool m_smoothing;
349 RealScalar m_angle;
350 bool m_residual;
351
352 public:
354 IDRS(): m_S(4), m_smoothing(false), m_angle(RealScalar(0.7)), m_residual(false) {}
355
366 template <typename MatrixDerived>
367 explicit IDRS(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_S(4), m_smoothing(false),
368 m_angle(RealScalar(0.7)), m_residual(false) {}
369
370
376 template <typename Rhs, typename Dest>
377 void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
378 {
379 m_iterations = Base::maxIterations();
380 m_error = Base::m_tolerance;
381
382 bool ret = internal::idrs(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error, m_S,m_smoothing,m_angle,m_residual);
383
384 m_info = (!ret) ? NumericalIssue : m_error <= Base::m_tolerance ? Success : NoConvergence;
385 }
386
388 void setS(Index S)
389 {
390 if (S < 1)
391 {
392 S = 4;
393 }
394
395 m_S = S;
396 }
397
404 void setSmoothing(bool smoothing)
405 {
406 m_smoothing=smoothing;
407 }
408
419 void setAngle(RealScalar angle)
420 {
421 m_angle=angle;
422 }
423
427 void setResidualUpdate(bool update)
428 {
429 m_residual=update;
430 }
431
432 };
433
434} // namespace Eigen
435
436#endif /* EIGEN_IDRS_H */
The Induced Dimension Reduction method (IDR(s)) is a short-recurrences Krylov method for sparse squar...
Definition: IDRS.h:332
void setSmoothing(bool smoothing)
Definition: IDRS.h:404
IDRS(const EigenBase< MatrixDerived > &A)
Definition: IDRS.h:367
void setS(Index S)
Definition: IDRS.h:388
void setResidualUpdate(bool update)
Definition: IDRS.h:427
IDRS()
Definition: IDRS.h:354
void _solve_vector_with_guess_impl(const Rhs &b, Dest &x) const
Definition: IDRS.h:377
void setAngle(RealScalar angle)
Definition: IDRS.h:419
NumericalIssue
Matrix< Type, Size, 1 > Vector
Namespace containing all symbols from the Eigen library.
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index