WPILibC++ 2023.4.3
BasicPreconditioners.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) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_BASIC_PRECONDITIONERS_H
11#define EIGEN_BASIC_PRECONDITIONERS_H
12
13namespace Eigen {
14
15/** \ingroup IterativeLinearSolvers_Module
16 * \brief A preconditioner based on the digonal entries
17 *
18 * This class allows to approximately solve for A.x = b problems assuming A is a diagonal matrix.
19 * In other words, this preconditioner neglects all off diagonal entries and, in Eigen's language, solves for:
20 \code
21 A.diagonal().asDiagonal() . x = b
22 \endcode
23 *
24 * \tparam _Scalar the type of the scalar.
25 *
26 * \implsparsesolverconcept
27 *
28 * This preconditioner is suitable for both selfadjoint and general problems.
29 * The diagonal entries are pre-inverted and stored into a dense vector.
30 *
31 * \note A variant that has yet to be implemented would attempt to preserve the norm of each column.
32 *
33 * \sa class LeastSquareDiagonalPreconditioner, class ConjugateGradient
34 */
35template <typename _Scalar>
37{
38 typedef _Scalar Scalar;
40 public:
41 typedef typename Vector::StorageIndex StorageIndex;
42 enum {
45 };
46
48
49 template<typename MatType>
50 explicit DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols())
51 {
52 compute(mat);
53 }
54
57
58 template<typename MatType>
60 {
61 return *this;
62 }
63
64 template<typename MatType>
65 DiagonalPreconditioner& factorize(const MatType& mat)
66 {
67 m_invdiag.resize(mat.cols());
68 for(int j=0; j<mat.outerSize(); ++j)
69 {
70 typename MatType::InnerIterator it(mat,j);
71 while(it && it.index()!=j) ++it;
72 if(it && it.index()==j && it.value()!=Scalar(0))
73 m_invdiag(j) = Scalar(1)/it.value();
74 else
75 m_invdiag(j) = Scalar(1);
76 }
77 m_isInitialized = true;
78 return *this;
79 }
80
81 template<typename MatType>
82 DiagonalPreconditioner& compute(const MatType& mat)
83 {
84 return factorize(mat);
85 }
86
87 /** \internal */
88 template<typename Rhs, typename Dest>
89 void _solve_impl(const Rhs& b, Dest& x) const
90 {
91 x = m_invdiag.array() * b.array() ;
92 }
93
94 template<typename Rhs> inline const Solve<DiagonalPreconditioner, Rhs>
95 solve(const MatrixBase<Rhs>& b) const
96 {
97 eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized.");
98 eigen_assert(m_invdiag.size()==b.rows()
99 && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b");
100 return Solve<DiagonalPreconditioner, Rhs>(*this, b.derived());
101 }
102
104
105 protected:
108};
109
110/** \ingroup IterativeLinearSolvers_Module
111 * \brief Jacobi preconditioner for LeastSquaresConjugateGradient
112 *
113 * This class allows to approximately solve for A' A x = A' b problems assuming A' A is a diagonal matrix.
114 * In other words, this preconditioner neglects all off diagonal entries and, in Eigen's language, solves for:
115 \code
116 (A.adjoint() * A).diagonal().asDiagonal() * x = b
117 \endcode
118 *
119 * \tparam _Scalar the type of the scalar.
120 *
121 * \implsparsesolverconcept
122 *
123 * The diagonal entries are pre-inverted and stored into a dense vector.
124 *
125 * \sa class LeastSquaresConjugateGradient, class DiagonalPreconditioner
126 */
127template <typename _Scalar>
129{
130 typedef _Scalar Scalar;
131 typedef typename NumTraits<Scalar>::Real RealScalar;
133 using Base::m_invdiag;
134 public:
135
137
138 template<typename MatType>
139 explicit LeastSquareDiagonalPreconditioner(const MatType& mat) : Base()
140 {
141 compute(mat);
142 }
143
144 template<typename MatType>
146 {
147 return *this;
148 }
149
150 template<typename MatType>
152 {
153 // Compute the inverse squared-norm of each column of mat
154 m_invdiag.resize(mat.cols());
155 if(MatType::IsRowMajor)
156 {
157 m_invdiag.setZero();
158 for(Index j=0; j<mat.outerSize(); ++j)
159 {
160 for(typename MatType::InnerIterator it(mat,j); it; ++it)
161 m_invdiag(it.index()) += numext::abs2(it.value());
162 }
163 for(Index j=0; j<mat.cols(); ++j)
164 if(numext::real(m_invdiag(j))>RealScalar(0))
165 m_invdiag(j) = RealScalar(1)/numext::real(m_invdiag(j));
166 }
167 else
168 {
169 for(Index j=0; j<mat.outerSize(); ++j)
170 {
171 RealScalar sum = mat.col(j).squaredNorm();
172 if(sum>RealScalar(0))
173 m_invdiag(j) = RealScalar(1)/sum;
174 else
175 m_invdiag(j) = RealScalar(1);
176 }
177 }
179 return *this;
180 }
181
182 template<typename MatType>
184 {
185 return factorize(mat);
186 }
187
189
190 protected:
191};
192
193/** \ingroup IterativeLinearSolvers_Module
194 * \brief A naive preconditioner which approximates any matrix as the identity matrix
195 *
196 * \implsparsesolverconcept
197 *
198 * \sa class DiagonalPreconditioner
199 */
201{
202 public:
203
205
206 template<typename MatrixType>
207 explicit IdentityPreconditioner(const MatrixType& ) {}
208
209 template<typename MatrixType>
210 IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; }
211
212 template<typename MatrixType>
213 IdentityPreconditioner& factorize(const MatrixType& ) { return *this; }
214
215 template<typename MatrixType>
216 IdentityPreconditioner& compute(const MatrixType& ) { return *this; }
217
218 template<typename Rhs>
219 inline const Rhs& solve(const Rhs& b) const { return b; }
220
222};
223
224} // end namespace Eigen
225
226#endif // EIGEN_BASIC_PRECONDITIONERS_H
EIGEN_DEVICE_FUNC RealReturnType real() const
Definition: CommonCwiseUnaryOps.h:100
#define EIGEN_NOEXCEPT
Definition: Macros.h:1428
#define EIGEN_CONSTEXPR
Definition: Macros.h:797
#define eigen_assert(x)
Definition: Macros.h:1047
A preconditioner based on the digonal entries.
Definition: BasicPreconditioners.h:37
bool m_isInitialized
Definition: BasicPreconditioners.h:107
Vector::StorageIndex StorageIndex
Definition: BasicPreconditioners.h:41
void _solve_impl(const Rhs &b, Dest &x) const
Definition: BasicPreconditioners.h:89
DiagonalPreconditioner & compute(const MatType &mat)
Definition: BasicPreconditioners.h:82
DiagonalPreconditioner & analyzePattern(const MatType &)
Definition: BasicPreconditioners.h:59
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: BasicPreconditioners.h:55
ComputationInfo info()
Definition: BasicPreconditioners.h:103
const Solve< DiagonalPreconditioner, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: BasicPreconditioners.h:95
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: BasicPreconditioners.h:56
Vector m_invdiag
Definition: BasicPreconditioners.h:106
@ MaxColsAtCompileTime
Definition: BasicPreconditioners.h:44
@ ColsAtCompileTime
Definition: BasicPreconditioners.h:43
DiagonalPreconditioner & factorize(const MatType &mat)
Definition: BasicPreconditioners.h:65
DiagonalPreconditioner(const MatType &mat)
Definition: BasicPreconditioners.h:50
DiagonalPreconditioner()
Definition: BasicPreconditioners.h:47
A naive preconditioner which approximates any matrix as the identity matrix.
Definition: BasicPreconditioners.h:201
IdentityPreconditioner & analyzePattern(const MatrixType &)
Definition: BasicPreconditioners.h:210
ComputationInfo info()
Definition: BasicPreconditioners.h:221
IdentityPreconditioner(const MatrixType &)
Definition: BasicPreconditioners.h:207
const Rhs & solve(const Rhs &b) const
Definition: BasicPreconditioners.h:219
IdentityPreconditioner & compute(const MatrixType &)
Definition: BasicPreconditioners.h:216
IdentityPreconditioner & factorize(const MatrixType &)
Definition: BasicPreconditioners.h:213
IdentityPreconditioner()
Definition: BasicPreconditioners.h:204
Jacobi preconditioner for LeastSquaresConjugateGradient.
Definition: BasicPreconditioners.h:129
LeastSquareDiagonalPreconditioner & factorize(const MatType &mat)
Definition: BasicPreconditioners.h:151
LeastSquareDiagonalPreconditioner()
Definition: BasicPreconditioners.h:136
ComputationInfo info()
Definition: BasicPreconditioners.h:188
LeastSquareDiagonalPreconditioner & analyzePattern(const MatType &)
Definition: BasicPreconditioners.h:145
LeastSquareDiagonalPreconditioner(const MatType &mat)
Definition: BasicPreconditioners.h:139
LeastSquareDiagonalPreconditioner & compute(const MatType &mat)
Definition: BasicPreconditioners.h:183
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Resizes *this to a rows x cols matrix.
Definition: PlainObjectBase.h:271
EIGEN_DEVICE_FUNC Derived & setZero(Index size)
Resizes to the given size, and sets all coefficients in this expression to zero.
Definition: CwiseNullaryOp.h:562
Pseudo expression representing a solving operation.
Definition: Solve.h:63
ComputationInfo
Enum for reporting the status of a computation.
Definition: Constants.h:440
@ Success
Computation was successful.
Definition: Constants.h:442
EIGEN_DEVICE_FUNC bool abs2(bool x)
Definition: MathFunctions.h:1292
Namespace containing all symbols from the Eigen library.
Definition: MatrixExponential.h:16
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const int Dynamic
This value means that a positive quantity (e.g., a size) is not known at compile-time,...
Definition: Constants.h:22
b
Definition: data.h:44
Holds information about the various numeric (i.e.
Definition: NumTraits.h:233