Eigen  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
BasicPreconditioners.h
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
13#include "./InternalHeaderCheck.h"
14
15namespace Eigen {
16
37template <typename Scalar_>
39{
40 typedef Scalar_ Scalar;
42 public:
43 typedef typename Vector::StorageIndex StorageIndex;
44 enum {
45 ColsAtCompileTime = Dynamic,
46 MaxColsAtCompileTime = Dynamic
47 };
48
49 DiagonalPreconditioner() : m_isInitialized(false) {}
50
51 template<typename MatType>
52 explicit DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols())
53 {
54 compute(mat);
55 }
56
57 EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_invdiag.size(); }
58 EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_invdiag.size(); }
59
60 template<typename MatType>
61 DiagonalPreconditioner& analyzePattern(const MatType& )
62 {
63 return *this;
64 }
65
66 template<typename MatType>
67 DiagonalPreconditioner& factorize(const MatType& mat)
68 {
69 m_invdiag.resize(mat.cols());
70 for(int j=0; j<mat.outerSize(); ++j)
71 {
72 typename MatType::InnerIterator it(mat,j);
73 while(it && it.index()!=j) ++it;
74 if(it && it.index()==j && it.value()!=Scalar(0))
75 m_invdiag(j) = Scalar(1)/it.value();
76 else
77 m_invdiag(j) = Scalar(1);
78 }
79 m_isInitialized = true;
80 return *this;
81 }
82
83 template<typename MatType>
84 DiagonalPreconditioner& compute(const MatType& mat)
85 {
86 return factorize(mat);
87 }
88
90 template<typename Rhs, typename Dest>
91 void _solve_impl(const Rhs& b, Dest& x) const
92 {
93 x = m_invdiag.array() * b.array() ;
94 }
95
96 template<typename Rhs> inline const Solve<DiagonalPreconditioner, Rhs>
97 solve(const MatrixBase<Rhs>& b) const
98 {
99 eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized.");
100 eigen_assert(m_invdiag.size()==b.rows()
101 && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b");
103 }
104
105 ComputationInfo info() { return Success; }
106
107 protected:
108 Vector m_invdiag;
109 bool m_isInitialized;
110};
111
129template <typename Scalar_>
131{
132 typedef Scalar_ Scalar;
133 typedef typename NumTraits<Scalar>::Real RealScalar;
135 using Base::m_invdiag;
136 public:
137
139
140 template<typename MatType>
141 explicit LeastSquareDiagonalPreconditioner(const MatType& mat) : Base()
142 {
143 compute(mat);
144 }
145
146 template<typename MatType>
147 LeastSquareDiagonalPreconditioner& analyzePattern(const MatType& )
148 {
149 return *this;
150 }
151
152 template<typename MatType>
153 LeastSquareDiagonalPreconditioner& factorize(const MatType& mat)
154 {
155 // Compute the inverse squared-norm of each column of mat
156 m_invdiag.resize(mat.cols());
157 if(MatType::IsRowMajor)
158 {
159 m_invdiag.setZero();
160 for(Index j=0; j<mat.outerSize(); ++j)
161 {
162 for(typename MatType::InnerIterator it(mat,j); it; ++it)
163 m_invdiag(it.index()) += numext::abs2(it.value());
164 }
165 for(Index j=0; j<mat.cols(); ++j)
166 if(numext::real(m_invdiag(j))>RealScalar(0))
167 m_invdiag(j) = RealScalar(1)/numext::real(m_invdiag(j));
168 }
169 else
170 {
171 for(Index j=0; j<mat.outerSize(); ++j)
172 {
173 RealScalar sum = mat.col(j).squaredNorm();
174 if(sum>RealScalar(0))
175 m_invdiag(j) = RealScalar(1)/sum;
176 else
177 m_invdiag(j) = RealScalar(1);
178 }
179 }
180 Base::m_isInitialized = true;
181 return *this;
182 }
183
184 template<typename MatType>
185 LeastSquareDiagonalPreconditioner& compute(const MatType& mat)
186 {
187 return factorize(mat);
188 }
189
190 ComputationInfo info() { return Success; }
191
192 protected:
193};
194
203{
204 public:
205
207
208 template<typename MatrixType>
209 explicit IdentityPreconditioner(const MatrixType& ) {}
210
211 template<typename MatrixType>
212 IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; }
213
214 template<typename MatrixType>
215 IdentityPreconditioner& factorize(const MatrixType& ) { return *this; }
216
217 template<typename MatrixType>
218 IdentityPreconditioner& compute(const MatrixType& ) { return *this; }
219
220 template<typename Rhs>
221 inline const Rhs& solve(const Rhs& b) const { return b; }
222
223 ComputationInfo info() { return Success; }
224};
225
226} // end namespace Eigen
227
228#endif // EIGEN_BASIC_PRECONDITIONERS_H
internal::traits< Derived >::StorageIndex StorageIndex
The type used to store indices.
Definition: DenseBase.h:58
Derived & derived()
Definition: EigenBase.h:48
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:62
A preconditioner based on the digonal entries.
Definition: BasicPreconditioners.h:39
A naive preconditioner which approximates any matrix as the identity matrix.
Definition: BasicPreconditioners.h:203
Jacobi preconditioner for LeastSquaresConjugateGradient.
Definition: BasicPreconditioners.h:131
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:283
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:564
Pseudo expression representing a solving operation.
Definition: Solve.h:65
ComputationInfo
Definition: Constants.h:442
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: B01_Experimental.dox:1
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const int Dynamic
Definition: Constants.h:24
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:235