Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
 
Loading...
Searching...
No Matches
IncompleteLU.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2011 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_INCOMPLETE_LU_H
11#define EIGEN_INCOMPLETE_LU_H
12
13namespace Eigen {
14
15template <typename _Scalar>
16class IncompleteLU : public SparseSolverBase<IncompleteLU<_Scalar> >
17{
18 protected:
19 typedef SparseSolverBase<IncompleteLU<_Scalar> > Base;
20 using Base::m_isInitialized;
21
22 typedef _Scalar Scalar;
23 typedef Matrix<Scalar,Dynamic,1> Vector;
24 typedef typename Vector::Index Index;
25 typedef SparseMatrix<Scalar,RowMajor> FactorType;
26
27 public:
28 typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
29
30 IncompleteLU() {}
31
32 template<typename MatrixType>
33 IncompleteLU(const MatrixType& mat)
34 {
35 compute(mat);
36 }
37
38 Index rows() const { return m_lu.rows(); }
39 Index cols() const { return m_lu.cols(); }
40
41 template<typename MatrixType>
42 IncompleteLU& compute(const MatrixType& mat)
43 {
44 m_lu = mat;
45 int size = mat.cols();
46 Vector diag(size);
47 for(int i=0; i<size; ++i)
48 {
49 typename FactorType::InnerIterator k_it(m_lu,i);
50 for(; k_it && k_it.index()<i; ++k_it)
51 {
52 int k = k_it.index();
53 k_it.valueRef() /= diag(k);
54
55 typename FactorType::InnerIterator j_it(k_it);
56 typename FactorType::InnerIterator kj_it(m_lu, k);
57 while(kj_it && kj_it.index()<=k) ++kj_it;
58 for(++j_it; j_it; )
59 {
60 if(kj_it.index()==j_it.index())
61 {
62 j_it.valueRef() -= k_it.value() * kj_it.value();
63 ++j_it;
64 ++kj_it;
65 }
66 else if(kj_it.index()<j_it.index()) ++kj_it;
67 else ++j_it;
68 }
69 }
70 if(k_it && k_it.index()==i) diag(i) = k_it.value();
71 else diag(i) = 1;
72 }
73 m_isInitialized = true;
74 return *this;
75 }
76
77 template<typename Rhs, typename Dest>
78 void _solve_impl(const Rhs& b, Dest& x) const
79 {
80 x = m_lu.template triangularView<UnitLower>().solve(b);
81 x = m_lu.template triangularView<Upper>().solve(x);
82 }
83
84 protected:
85 FactorType m_lu;
86};
87
88} // end namespace Eigen
89
90#endif // EIGEN_INCOMPLETE_LU_H
Matrix< Type, Size, 1 > Vector
Namespace containing all symbols from the Eigen library.
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index