IDRS.h
Go to the documentation of this file.
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 
16 namespace 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;
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 
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;
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;
349  RealScalar m_angle;
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 
428  {
429  m_residual=update;
430  }
431 
432  };
433 
434 } // namespace Eigen
435 
436 #endif /* EIGEN_IDRS_H */
bool m_residual
Definition: IDRS.h:350
const Solve< FullPivLU< _MatrixType >, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SolverBase.h:106
SCALAR Scalar
Definition: bench_gemm.cpp:46
MatrixType::RealScalar RealScalar
Definition: IDRS.h:337
Matrix< RealScalar, Dynamic, Dynamic > M
Definition: bench_gemm.cpp:51
void setSmoothing(bool smoothing)
Definition: IDRS.h:404
Scalar * b
Definition: benchVecAdd.cpp:17
void adjoint(const MatrixType &m)
Definition: adjoint.cpp:67
IDRS()
Definition: IDRS.h:354
def update(text)
Definition: relicense.py:46
JacobiRotation< float > G
FullPivLU & compute(const EigenBase< InputType > &matrix)
Definition: FullPivLU.h:120
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
MatrixXf MatrixType
iterator iter(handle obj)
Definition: pytypes.h:2273
The Induced Dimension Reduction method (IDR(s)) is a short-recurrences Krylov method for sparse squar...
Definition: IDRS.h:275
#define N
Definition: gksort.c:12
DiscreteKey S(1, 2)
HouseholderQR< MatrixXf > qr(A)
static double epsilon
Definition: testRot3.cpp:37
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE internal::enable_if< NumTraits< T >::IsSigned||NumTraits< T >::IsComplex, typename NumTraits< T >::Real >::type abs(const T &x)
bool m_smoothing
Definition: IDRS.h:348
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Array< int, Dynamic, 1 > v
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
RealScalar alpha
Vector::Scalar omega(const Vector &t, const Vector &s, RealScalar angle)
Definition: IDRS.h:36
RealScalar s
IterativeSolverBase< IDRS > Base
Definition: IDRS.h:341
HouseholderSequenceType householderQ() const
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
DenseIndex ret
void _solve_vector_with_guess_impl(const Rhs &b, Dest &x) const
Definition: IDRS.h:377
MatrixType::Scalar Scalar
Definition: IDRS.h:336
_MatrixType MatrixType
Definition: IDRS.h:335
LU decomposition of a matrix with complete pivoting, and related features.
void setResidualUpdate(bool update)
Definition: IDRS.h:427
Householder QR decomposition of a matrix.
RealScalar m_angle
Definition: IDRS.h:349
void setS(Index S)
Definition: IDRS.h:388
Index m_S
Definition: IDRS.h:347
IDRS(const EigenBase< MatrixDerived > &A)
Definition: IDRS.h:367
const G double tol
Definition: Group.h:86
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
EIGEN_DEVICE_FUNC bool isApprox(const Scalar &x, const Scalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
#define abs(x)
Definition: datatypes.h:17
Eigen::Matrix< double, Eigen::Dynamic, 1 > Vector
Base class for linear iterative solvers.
void setAngle(RealScalar angle)
Definition: IDRS.h:419
_Preconditioner Preconditioner
Definition: IDRS.h:338
bool idrs(const MatrixType &A, const Rhs &b, Dest &x, const Preconditioner &precond, Index &iter, typename Dest::RealScalar &relres, Index S, bool smoothing, typename Dest::RealScalar angle, bool replacement)
Definition: IDRS.h:59
Point2 t(10, 10)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:21