MINRES.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) 2012 Giacomo Po <gpo@ucla.edu>
5 // Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
6 // Copyright (C) 2018 David Hyde <dabh@stanford.edu>
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_MINRES_H_
14 #define EIGEN_MINRES_H_
15 
16 
17 namespace Eigen {
18 
19  namespace internal {
20 
30  template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
32  void minres(const MatrixType& mat, const Rhs& rhs, Dest& x,
33  const Preconditioner& precond, Index& iters,
34  typename Dest::RealScalar& tol_error)
35  {
36  using std::sqrt;
37  typedef typename Dest::RealScalar RealScalar;
38  typedef typename Dest::Scalar Scalar;
40 
41  // Check for zero rhs
42  const RealScalar rhsNorm2(rhs.squaredNorm());
43  if(rhsNorm2 == 0)
44  {
45  x.setZero();
46  iters = 0;
47  tol_error = 0;
48  return;
49  }
50 
51  // initialize
52  const Index maxIters(iters); // initialize maxIters to iters
53  const Index N(mat.cols()); // the size of the matrix
54  const RealScalar threshold2(tol_error*tol_error*rhsNorm2); // convergence threshold (compared to residualNorm2)
55 
56  // Initialize preconditioned Lanczos
57  VectorType v_old(N); // will be initialized inside loop
58  VectorType v( VectorType::Zero(N) ); //initialize v
59  VectorType v_new(rhs-mat*x); //initialize v_new
60  RealScalar residualNorm2(v_new.squaredNorm());
61  VectorType w(N); // will be initialized inside loop
62  VectorType w_new(precond.solve(v_new)); // initialize w_new
63 // RealScalar beta; // will be initialized inside loop
64  RealScalar beta_new2(v_new.dot(w_new));
65  eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
66  RealScalar beta_new(sqrt(beta_new2));
67  const RealScalar beta_one(beta_new);
68  // Initialize other variables
69  RealScalar c(1.0); // the cosine of the Givens rotation
70  RealScalar c_old(1.0);
71  RealScalar s(0.0); // the sine of the Givens rotation
72  RealScalar s_old(0.0); // the sine of the Givens rotation
73  VectorType p_oold(N); // will be initialized in loop
74  VectorType p_old(VectorType::Zero(N)); // initialize p_old=0
75  VectorType p(p_old); // initialize p=0
76  RealScalar eta(1.0);
77 
78  iters = 0; // reset iters
79  while ( iters < maxIters )
80  {
81  // Preconditioned Lanczos
82  /* Note that there are 4 variants on the Lanczos algorithm. These are
83  * described in Paige, C. C. (1972). Computational variants of
84  * the Lanczos method for the eigenproblem. IMA Journal of Applied
85  * Mathematics, 10(3), 373-381. The current implementation corresponds
86  * to the case A(2,7) in the paper. It also corresponds to
87  * algorithm 6.14 in Y. Saad, Iterative Methods for Sparse Linear
88  * Systems, 2003 p.173. For the preconditioned version see
89  * A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM (1987).
90  */
91  const RealScalar beta(beta_new);
92  v_old = v; // update: at first time step, this makes v_old = 0 so value of beta doesn't matter
93  v_new /= beta_new; // overwrite v_new for next iteration
94  w_new /= beta_new; // overwrite w_new for next iteration
95  v = v_new; // update
96  w = w_new; // update
97  v_new.noalias() = mat*w - beta*v_old; // compute v_new
98  const RealScalar alpha = v_new.dot(w);
99  v_new -= alpha*v; // overwrite v_new
100  w_new = precond.solve(v_new); // overwrite w_new
101  beta_new2 = v_new.dot(w_new); // compute beta_new
102  eigen_assert(beta_new2 >= 0.0 && "PRECONDITIONER IS NOT POSITIVE DEFINITE");
103  beta_new = sqrt(beta_new2); // compute beta_new
104 
105  // Givens rotation
106  const RealScalar r2 =s*alpha+c*c_old*beta; // s, s_old, c and c_old are still from previous iteration
107  const RealScalar r3 =s_old*beta; // s, s_old, c and c_old are still from previous iteration
108  const RealScalar r1_hat=c*alpha-c_old*s*beta;
109  const RealScalar r1 =sqrt( std::pow(r1_hat,2) + std::pow(beta_new,2) );
110  c_old = c; // store for next iteration
111  s_old = s; // store for next iteration
112  c=r1_hat/r1; // new cosine
113  s=beta_new/r1; // new sine
114 
115  // Update solution
116  p_oold = p_old;
117  p_old = p;
118  p.noalias()=(w-r2*p_old-r3*p_oold) /r1; // IS NOALIAS REQUIRED?
119  x += beta_one*c*eta*p;
120 
121  /* Update the squared residual. Note that this is the estimated residual.
122  The real residual |Ax-b|^2 may be slightly larger */
123  residualNorm2 *= s*s;
124 
125  if ( residualNorm2 < threshold2)
126  {
127  break;
128  }
129 
130  eta=-s*eta; // update eta
131  iters++; // increment iteration number (for output purposes)
132  }
133 
134  /* Compute error. Note that this is the estimated error. The real
135  error |Ax-b|/|b| may be slightly larger */
136  tol_error = std::sqrt(residualNorm2 / rhsNorm2);
137  }
138 
139  }
140 
141  template< typename _MatrixType, int _UpLo=Lower,
142  typename _Preconditioner = IdentityPreconditioner>
143  class MINRES;
144 
145  namespace internal {
146 
147  template< typename _MatrixType, int _UpLo, typename _Preconditioner>
148  struct traits<MINRES<_MatrixType,_UpLo,_Preconditioner> >
149  {
150  typedef _MatrixType MatrixType;
151  typedef _Preconditioner Preconditioner;
152  };
153 
154  }
155 
194  template< typename _MatrixType, int _UpLo, typename _Preconditioner>
195  class MINRES : public IterativeSolverBase<MINRES<_MatrixType,_UpLo,_Preconditioner> >
196  {
197 
199  using Base::matrix;
200  using Base::m_error;
201  using Base::m_iterations;
202  using Base::m_info;
203  using Base::m_isInitialized;
204  public:
205  using Base::_solve_impl;
206  typedef _MatrixType MatrixType;
207  typedef typename MatrixType::Scalar Scalar;
209  typedef _Preconditioner Preconditioner;
210 
211  enum {UpLo = _UpLo};
212 
213  public:
214 
216  MINRES() : Base() {}
217 
228  template<typename MatrixDerived>
229  explicit MINRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
230 
233 
235  template<typename Rhs,typename Dest>
236  void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
237  {
238  typedef typename Base::MatrixWrapper MatrixWrapper;
239  typedef typename Base::ActualMatrixType ActualMatrixType;
240  enum {
241  TransposeInput = (!MatrixWrapper::MatrixFree)
242  && (UpLo==(Lower|Upper))
243  && (!MatrixType::IsRowMajor)
245  };
246  typedef typename internal::conditional<TransposeInput,Transpose<const ActualMatrixType>, ActualMatrixType const&>::type RowMajorWrapper;
247  EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(MatrixWrapper::MatrixFree,UpLo==(Lower|Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY);
248  typedef typename internal::conditional<UpLo==(Lower|Upper),
249  RowMajorWrapper,
250  typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
251  >::type SelfAdjointWrapper;
252 
255  RowMajorWrapper row_mat(matrix());
256  internal::minres(SelfAdjointWrapper(row_mat), b, x,
259  }
260 
261  protected:
262 
263  };
264 
265 } // end namespace Eigen
266 
267 #endif // EIGEN_MINRES_H
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Eigen::MINRES::MINRES
MINRES(const EigenBase< MatrixDerived > &A)
Definition: MINRES.h:229
Eigen::IterativeSolverBase::_solve_impl
void _solve_impl(const Rhs &b, Dest &x) const
Definition: IterativeSolverBase.h:400
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:147
s
RealScalar s
Definition: level1_cplx_impl.h:126
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
r2
static const double r2
Definition: testSmartRangeFactor.cpp:32
Eigen::internal::traits< MINRES< _MatrixType, _UpLo, _Preconditioner > >::MatrixType
_MatrixType MatrixType
Definition: MINRES.h:150
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
b
Scalar * b
Definition: benchVecAdd.cpp:17
Eigen::EigenBase
Definition: EigenBase.h:29
eigen_assert
#define eigen_assert(x)
Definition: Macros.h:1037
x
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
Definition: gnuplot_common_settings.hh:12
Eigen::Upper
@ Upper
Definition: Constants.h:211
Eigen::Success
@ Success
Definition: Constants.h:442
type
Definition: pytypes.h:1525
mat
MatrixXf mat
Definition: Tutorial_AdvancedInitialization_CommaTemporary.cpp:1
Eigen::MINRES::UpLo
@ UpLo
Definition: MINRES.h:211
r1
static const double r1
Definition: testSmartRangeFactor.cpp:32
beta
double beta(double a, double b)
Definition: beta.c:61
Eigen::NoConvergence
@ NoConvergence
Definition: Constants.h:446
A
Definition: test_numpy_dtypes.cpp:298
Eigen::MINRES
A minimal residual solver for sparse symmetric problems.
Definition: MINRES.h:143
Eigen::MINRES::Preconditioner
_Preconditioner Preconditioner
Definition: MINRES.h:209
Eigen::internal::traits< MINRES< _MatrixType, _UpLo, _Preconditioner > >::Preconditioner
_Preconditioner Preconditioner
Definition: MINRES.h:151
Eigen::Architecture::Type
Type
Definition: Constants.h:471
Eigen::IterativeSolverBase::maxIterations
Index maxIterations() const
Definition: IterativeSolverBase.h:281
r3
static const double r3
Definition: testSmartRangeFactor.cpp:32
Eigen::IterativeSolverBase::m_info
ComputationInfo m_info
Definition: IterativeSolverBase.h:438
ceres::pow
Jet< T, N > pow(const Jet< T, N > &f, double g)
Definition: jet.h:570
Eigen::Lower
@ Lower
Definition: Constants.h:209
matrix
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Definition: gtsam/3rdparty/Eigen/blas/common.h:110
Eigen::IterativeSolverBase
Base class for linear iterative solvers.
Definition: IterativeSolverBase.h:143
Eigen::MINRES::m_iterations
Index m_iterations
Definition: IterativeSolverBase.h:437
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
Eigen::MINRES::_solve_vector_with_guess_impl
void _solve_vector_with_guess_impl(const Rhs &b, Dest &x) const
Definition: MINRES.h:236
Eigen::IterativeSolverBase::m_tolerance
RealScalar m_tolerance
Definition: IterativeSolverBase.h:434
Eigen::internal::traits
Definition: ForwardDeclarations.h:17
EIGEN_STATIC_ASSERT
#define EIGEN_STATIC_ASSERT(CONDITION, MSG)
Definition: StaticAssert.h:127
Eigen::MINRES::MINRES
MINRES()
Definition: MINRES.h:216
Eigen::internal::Rhs
@ Rhs
Definition: TensorContractionMapper.h:18
Eigen::MINRES::~MINRES
~MINRES()
Definition: MINRES.h:232
Eigen::MINRES::m_error
RealScalar m_error
Definition: IterativeSolverBase.h:436
Eigen::internal::conditional
Definition: Meta.h:109
Eigen::IterativeSolverBase::m_iterations
Index m_iterations
Definition: IterativeSolverBase.h:437
p
float * p
Definition: Tutorial_Map_using.cpp:9
Eigen::IterativeSolverBase::matrix
const ActualMatrixType & matrix() const
Definition: IterativeSolverBase.h:419
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
Eigen::internal::minres
EIGEN_DONT_INLINE void minres(const MatrixType &mat, const Rhs &rhs, Dest &x, const Preconditioner &precond, Index &iters, typename Dest::RealScalar &tol_error)
Definition: MINRES.h:32
Eigen::Matrix< Scalar, Dynamic, 1 >
Eigen::MINRES::RealScalar
MatrixType::RealScalar RealScalar
Definition: MINRES.h:208
Eigen::MINRES::Base
IterativeSolverBase< MINRES > Base
Definition: MINRES.h:198
N
#define N
Definition: igam.h:9
internal
Definition: BandTriangularSolver.h:13
VectorType
Definition: FFTW.cpp:65
EIGEN_IMPLIES
#define EIGEN_IMPLIES(a, b)
Definition: Macros.h:1315
Eigen::MatrixWrapper
Expression of an array as a mathematical vector or matrix.
Definition: ArrayBase.h:15
Eigen::internal::generic_matrix_wrapper< MatrixType >
Eigen::IterativeSolverBase::ActualMatrixType
MatrixWrapper::ActualMatrixType ActualMatrixType
Definition: IterativeSolverBase.h:417
Eigen::IterativeSolverBase::m_preconditioner
Preconditioner m_preconditioner
Definition: IterativeSolverBase.h:431
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:232
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
Eigen::MINRES::Scalar
MatrixType::Scalar Scalar
Definition: MINRES.h:207
Eigen::MINRES::m_info
ComputationInfo m_info
Definition: IterativeSolverBase.h:438
Eigen::IterativeSolverBase::m_error
RealScalar m_error
Definition: IterativeSolverBase.h:436
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
Eigen::MINRES::MatrixType
_MatrixType MatrixType
Definition: MINRES.h:206
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
EIGEN_DONT_INLINE
#define EIGEN_DONT_INLINE
Definition: Macros.h:940
Eigen::SparseSolverBase::m_isInitialized
bool m_isInitialized
Definition: SparseSolverBase.h:119


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:33:49