Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
KLUSupport.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2017 Kyle Macfarlan <kyle.macfarlan@gmail.com>
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_KLUSUPPORT_H
11#define EIGEN_KLUSUPPORT_H
12
13namespace Eigen {
14
15/* TODO extract L, extract U, compute det, etc... */
16
34inline int klu_solve(klu_symbolic *Symbolic, klu_numeric *Numeric, Index ldim, Index nrhs, double B [ ], klu_common *Common, double) {
35 return klu_solve(Symbolic, Numeric, internal::convert_index<int>(ldim), internal::convert_index<int>(nrhs), B, Common);
36}
37
38inline int klu_solve(klu_symbolic *Symbolic, klu_numeric *Numeric, Index ldim, Index nrhs, std::complex<double>B[], klu_common *Common, std::complex<double>) {
39 return klu_z_solve(Symbolic, Numeric, internal::convert_index<int>(ldim), internal::convert_index<int>(nrhs), &numext::real_ref(B[0]), Common);
40}
41
42inline int klu_tsolve(klu_symbolic *Symbolic, klu_numeric *Numeric, Index ldim, Index nrhs, double B[], klu_common *Common, double) {
43 return klu_tsolve(Symbolic, Numeric, internal::convert_index<int>(ldim), internal::convert_index<int>(nrhs), B, Common);
44}
45
46inline int klu_tsolve(klu_symbolic *Symbolic, klu_numeric *Numeric, Index ldim, Index nrhs, std::complex<double>B[], klu_common *Common, std::complex<double>) {
47 return klu_z_tsolve(Symbolic, Numeric, internal::convert_index<int>(ldim), internal::convert_index<int>(nrhs), &numext::real_ref(B[0]), 0, Common);
48}
49
50inline klu_numeric* klu_factor(int Ap [ ], int Ai [ ], double Ax [ ], klu_symbolic *Symbolic, klu_common *Common, double) {
51 return klu_factor(Ap, Ai, Ax, Symbolic, Common);
52}
53
54inline klu_numeric* klu_factor(int Ap[], int Ai[], std::complex<double> Ax[], klu_symbolic *Symbolic, klu_common *Common, std::complex<double>) {
55 return klu_z_factor(Ap, Ai, &numext::real_ref(Ax[0]), Symbolic, Common);
56}
57
58
59template<typename _MatrixType>
60class KLU : public SparseSolverBase<KLU<_MatrixType> >
61{
62 protected:
63 typedef SparseSolverBase<KLU<_MatrixType> > Base;
64 using Base::m_isInitialized;
65 public:
66 using Base::_solve_impl;
67 typedef _MatrixType MatrixType;
68 typedef typename MatrixType::Scalar Scalar;
69 typedef typename MatrixType::RealScalar RealScalar;
70 typedef typename MatrixType::StorageIndex StorageIndex;
71 typedef Matrix<Scalar,Dynamic,1> Vector;
72 typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
73 typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
74 typedef SparseMatrix<Scalar> LUMatrixType;
75 typedef SparseMatrix<Scalar,ColMajor,int> KLUMatrixType;
76 typedef Ref<const KLUMatrixType, StandardCompressedFormat> KLUMatrixRef;
77 enum {
78 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
79 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
80 };
81
82 public:
83
84 KLU()
85 : m_dummy(0,0), mp_matrix(m_dummy)
86 {
87 init();
88 }
89
90 template<typename InputMatrixType>
91 explicit KLU(const InputMatrixType& matrix)
92 : mp_matrix(matrix)
93 {
94 init();
95 compute(matrix);
96 }
97
98 ~KLU()
99 {
100 if(m_symbolic) klu_free_symbolic(&m_symbolic,&m_common);
101 if(m_numeric) klu_free_numeric(&m_numeric,&m_common);
102 }
103
104 EIGEN_CONSTEXPR inline Index rows() const EIGEN_NOEXCEPT { return mp_matrix.rows(); }
105 EIGEN_CONSTEXPR inline Index cols() const EIGEN_NOEXCEPT { return mp_matrix.cols(); }
106
112 ComputationInfo info() const
113 {
114 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
115 return m_info;
116 }
117#if 0 // not implemented yet
118 inline const LUMatrixType& matrixL() const
119 {
120 if (m_extractedDataAreDirty) extractData();
121 return m_l;
122 }
123
124 inline const LUMatrixType& matrixU() const
125 {
126 if (m_extractedDataAreDirty) extractData();
127 return m_u;
128 }
129
130 inline const IntColVectorType& permutationP() const
131 {
132 if (m_extractedDataAreDirty) extractData();
133 return m_p;
134 }
135
136 inline const IntRowVectorType& permutationQ() const
137 {
138 if (m_extractedDataAreDirty) extractData();
139 return m_q;
140 }
141#endif
146 template<typename InputMatrixType>
147 void compute(const InputMatrixType& matrix)
148 {
149 if(m_symbolic) klu_free_symbolic(&m_symbolic, &m_common);
150 if(m_numeric) klu_free_numeric(&m_numeric, &m_common);
151 grab(matrix.derived());
152 analyzePattern_impl();
153 factorize_impl();
154 }
155
162 template<typename InputMatrixType>
163 void analyzePattern(const InputMatrixType& matrix)
164 {
165 if(m_symbolic) klu_free_symbolic(&m_symbolic, &m_common);
166 if(m_numeric) klu_free_numeric(&m_numeric, &m_common);
167
168 grab(matrix.derived());
169
170 analyzePattern_impl();
171 }
172
173
178 inline const klu_common& kluCommon() const
179 {
180 return m_common;
181 }
182
189 inline klu_common& kluCommon()
190 {
191 return m_common;
192 }
193
200 template<typename InputMatrixType>
201 void factorize(const InputMatrixType& matrix)
202 {
203 eigen_assert(m_analysisIsOk && "KLU: you must first call analyzePattern()");
204 if(m_numeric)
205 klu_free_numeric(&m_numeric,&m_common);
206
207 grab(matrix.derived());
208
209 factorize_impl();
210 }
211
213 template<typename BDerived,typename XDerived>
214 bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const;
215
216#if 0 // not implemented yet
217 Scalar determinant() const;
218
219 void extractData() const;
220#endif
221
222 protected:
223
224 void init()
225 {
226 m_info = InvalidInput;
227 m_isInitialized = false;
228 m_numeric = 0;
229 m_symbolic = 0;
230 m_extractedDataAreDirty = true;
231
232 klu_defaults(&m_common);
233 }
234
235 void analyzePattern_impl()
236 {
237 m_info = InvalidInput;
238 m_analysisIsOk = false;
239 m_factorizationIsOk = false;
240 m_symbolic = klu_analyze(internal::convert_index<int>(mp_matrix.rows()),
241 const_cast<StorageIndex*>(mp_matrix.outerIndexPtr()), const_cast<StorageIndex*>(mp_matrix.innerIndexPtr()),
242 &m_common);
243 if (m_symbolic) {
244 m_isInitialized = true;
245 m_info = Success;
246 m_analysisIsOk = true;
247 m_extractedDataAreDirty = true;
248 }
249 }
250
251 void factorize_impl()
252 {
253
254 m_numeric = klu_factor(const_cast<StorageIndex*>(mp_matrix.outerIndexPtr()), const_cast<StorageIndex*>(mp_matrix.innerIndexPtr()), const_cast<Scalar*>(mp_matrix.valuePtr()),
255 m_symbolic, &m_common, Scalar());
256
257
258 m_info = m_numeric ? Success : NumericalIssue;
259 m_factorizationIsOk = m_numeric ? 1 : 0;
260 m_extractedDataAreDirty = true;
261 }
262
263 template<typename MatrixDerived>
264 void grab(const EigenBase<MatrixDerived> &A)
265 {
266 mp_matrix.~KLUMatrixRef();
267 ::new (&mp_matrix) KLUMatrixRef(A.derived());
268 }
269
270 void grab(const KLUMatrixRef &A)
271 {
272 if(&(A.derived()) != &mp_matrix)
273 {
274 mp_matrix.~KLUMatrixRef();
275 ::new (&mp_matrix) KLUMatrixRef(A);
276 }
277 }
278
279 // cached data to reduce reallocation, etc.
280#if 0 // not implemented yet
281 mutable LUMatrixType m_l;
282 mutable LUMatrixType m_u;
283 mutable IntColVectorType m_p;
284 mutable IntRowVectorType m_q;
285#endif
286
287 KLUMatrixType m_dummy;
288 KLUMatrixRef mp_matrix;
289
290 klu_numeric* m_numeric;
291 klu_symbolic* m_symbolic;
292 klu_common m_common;
293 mutable ComputationInfo m_info;
294 int m_factorizationIsOk;
295 int m_analysisIsOk;
296 mutable bool m_extractedDataAreDirty;
297
298 private:
299 KLU(const KLU& ) { }
300};
301
302#if 0 // not implemented yet
303template<typename MatrixType>
304void KLU<MatrixType>::extractData() const
305{
306 if (m_extractedDataAreDirty)
307 {
308 eigen_assert(false && "KLU: extractData Not Yet Implemented");
309
310 // get size of the data
311 int lnz, unz, rows, cols, nz_udiag;
312 umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
313
314 // allocate data
315 m_l.resize(rows,(std::min)(rows,cols));
316 m_l.resizeNonZeros(lnz);
317
318 m_u.resize((std::min)(rows,cols),cols);
319 m_u.resizeNonZeros(unz);
320
321 m_p.resize(rows);
322 m_q.resize(cols);
323
324 // extract
325 umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(),
326 m_u.outerIndexPtr(), m_u.innerIndexPtr(), m_u.valuePtr(),
327 m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
328
329 m_extractedDataAreDirty = false;
330 }
331}
332
333template<typename MatrixType>
334typename KLU<MatrixType>::Scalar KLU<MatrixType>::determinant() const
335{
336 eigen_assert(false && "KLU: extractData Not Yet Implemented");
337 return Scalar();
338}
339#endif
340
341template<typename MatrixType>
342template<typename BDerived,typename XDerived>
343bool KLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const
344{
345 Index rhsCols = b.cols();
346 EIGEN_STATIC_ASSERT((XDerived::Flags&RowMajorBit)==0, THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
347 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
348
349 x = b;
350 int info = klu_solve(m_symbolic, m_numeric, b.rows(), rhsCols, x.const_cast_derived().data(), const_cast<klu_common*>(&m_common), Scalar());
351
352 m_info = info!=0 ? Success : NumericalIssue;
353 return true;
354}
355
356} // end namespace Eigen
357
358#endif // EIGEN_KLUSUPPORT_H
int klu_solve(klu_symbolic *Symbolic, klu_numeric *Numeric, Index ldim, Index nrhs, double B[], klu_common *Common, double)
A sparse LU factorization and solver based on KLU.
Definition: KLUSupport.h:34
ComputationInfo
Definition: Constants.h:440
@ NumericalIssue
Definition: Constants.h:444
@ InvalidInput
Definition: Constants.h:449
@ Success
Definition: Constants.h:442
Matrix< Type, Size, 1 > Vector
[c++11]
Definition: Matrix.h:551
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74