Context Navigation

source: branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Eigen/Eigen/src/SuperLUSupport/SuperLUSupport.h @ 9562

Visit:

Last change on this file since 9562 was 9562, checked in by gkronber, 11 years ago
#1967 worked on Gaussian process evolution.
File size: 31.8 KB

Line
1	// This file is part of Eigen, a lightweight C++ template library
2	// for linear algebra.
3	//
4	// Copyright (C) 2008-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_SUPERLUSUPPORT_H
11	#define EIGEN_SUPERLUSUPPORT_H
12
13	namespace Eigen {
14
15	#define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE) \
16	extern "C" { \
17	typedef struct { FLOATTYPE for_lu; FLOATTYPE total_needed; int expansions; } PREFIX##mem_usage_t; \
18	extern void PREFIX##gssvx(superlu_options_t , SuperMatrix , int , int , int *, \
19	char , FLOATTYPE , FLOATTYPE , SuperMatrix , SuperMatrix *, \
20	void , int, SuperMatrix , SuperMatrix *, \
21	FLOATTYPE , FLOATTYPE , FLOATTYPE , FLOATTYPE , \
22	PREFIX##mem_usage_t , SuperLUStat_t , int *); \
23	} \
24	inline float SuperLU_gssvx(superlu_options_t options, SuperMatrix A, \
25	int perm_c, int perm_r, int etree, char equed, \
26	FLOATTYPE R, FLOATTYPE C, SuperMatrix *L, \
27	SuperMatrix U, void work, int lwork, \
28	SuperMatrix B, SuperMatrix X, \
29	FLOATTYPE *recip_pivot_growth, \
30	FLOATTYPE rcond, FLOATTYPE ferr, FLOATTYPE *berr, \
31	SuperLUStat_t stats, int info, KEYTYPE) { \
32	PREFIX##mem_usage_t mem_usage; \
33	PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
34	U, work, lwork, B, X, recip_pivot_growth, rcond, \
35	ferr, berr, &mem_usage, stats, info); \
36	return mem_usage.for_lu; /* bytes used by the factor storage */ \
37	}
38
39	DECL_GSSVX(s,float,float)
40	DECL_GSSVX(c,float,std::complex<float>)
41	DECL_GSSVX(d,double,double)
42	DECL_GSSVX(z,double,std::complex<double>)
43
44	#ifdef MILU_ALPHA
45	#define EIGEN_SUPERLU_HAS_ILU
46	#endif
47
48	#ifdef EIGEN_SUPERLU_HAS_ILU
49
50	// similarly for the incomplete factorization using gsisx
51	#define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE) \
52	extern "C" { \
53	extern void PREFIX##gsisx(superlu_options_t , SuperMatrix , int , int , int *, \
54	char , FLOATTYPE , FLOATTYPE , SuperMatrix , SuperMatrix *, \
55	void , int, SuperMatrix , SuperMatrix , FLOATTYPE , FLOATTYPE *, \
56	PREFIX##mem_usage_t , SuperLUStat_t , int *); \
57	} \
58	inline float SuperLU_gsisx(superlu_options_t options, SuperMatrix A, \
59	int perm_c, int perm_r, int etree, char equed, \
60	FLOATTYPE R, FLOATTYPE C, SuperMatrix *L, \
61	SuperMatrix U, void work, int lwork, \
62	SuperMatrix B, SuperMatrix X, \
63	FLOATTYPE *recip_pivot_growth, \
64	FLOATTYPE *rcond, \
65	SuperLUStat_t stats, int info, KEYTYPE) { \
66	PREFIX##mem_usage_t mem_usage; \
67	PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
68	U, work, lwork, B, X, recip_pivot_growth, rcond, \
69	&mem_usage, stats, info); \
70	return mem_usage.for_lu; /* bytes used by the factor storage */ \
71	}
72
73	DECL_GSISX(s,float,float)
74	DECL_GSISX(c,float,std::complex<float>)
75	DECL_GSISX(d,double,double)
76	DECL_GSISX(z,double,std::complex<double>)
77
78	#endif
79
80	template<typename MatrixType>
81	struct SluMatrixMapHelper;
82
83	/** \internal
84	*
85	* A wrapper class for SuperLU matrices. It supports only compressed sparse matrices
86	* and dense matrices. Supernodal and other fancy format are not supported by this wrapper.
87	*
88	* This wrapper class mainly aims to avoids the need of dynamic allocation of the storage structure.
89	*/
90	struct SluMatrix : SuperMatrix
91	{
92	SluMatrix()
93	{
94	Store = &storage;
95	}
96
97	SluMatrix(const SluMatrix& other)
98	: SuperMatrix(other)
99	{
100	Store = &storage;
101	storage = other.storage;
102	}
103
104	SluMatrix& operator=(const SluMatrix& other)
105	{
106	SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
107	Store = &storage;
108	storage = other.storage;
109	return *this;
110	}
111
112	struct
113	{
114	union {int nnz;int lda;};
115	void *values;
116	int *innerInd;
117	int *outerInd;
118	} storage;
119
120	void setStorageType(Stype_t t)
121	{
122	Stype = t;
123	if (t==SLU_NC \|\| t==SLU_NR \|\| t==SLU_DN)
124	Store = &storage;
125	else
126	{
127	eigen_assert(false && "storage type not supported");
128	Store = 0;
129	}
130	}
131
132	template<typename Scalar>
133	void setScalarType()
134	{
135	if (internal::is_same<Scalar,float>::value)
136	Dtype = SLU_S;
137	else if (internal::is_same<Scalar,double>::value)
138	Dtype = SLU_D;
139	else if (internal::is_same<Scalar,std::complex<float> >::value)
140	Dtype = SLU_C;
141	else if (internal::is_same<Scalar,std::complex<double> >::value)
142	Dtype = SLU_Z;
143	else
144	{
145	eigen_assert(false && "Scalar type not supported by SuperLU");
146	}
147	}
148
149	template<typename MatrixType>
150	static SluMatrix Map(MatrixBase<MatrixType>& _mat)
151	{
152	MatrixType& mat(_mat.derived());
153	eigen_assert( ((MatrixType::Flags&RowMajorBit)!=RowMajorBit) && "row-major dense matrices are not supported by SuperLU");
154	SluMatrix res;
155	res.setStorageType(SLU_DN);
156	res.setScalarType<typename MatrixType::Scalar>();
157	res.Mtype = SLU_GE;
158
159	res.nrow = mat.rows();
160	res.ncol = mat.cols();
161
162	res.storage.lda = MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride();
163	res.storage.values = mat.data();
164	return res;
165	}
166
167	template<typename MatrixType>
168	static SluMatrix Map(SparseMatrixBase<MatrixType>& mat)
169	{
170	SluMatrix res;
171	if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
172	{
173	res.setStorageType(SLU_NR);
174	res.nrow = mat.cols();
175	res.ncol = mat.rows();
176	}
177	else
178	{
179	res.setStorageType(SLU_NC);
180	res.nrow = mat.rows();
181	res.ncol = mat.cols();
182	}
183
184	res.Mtype = SLU_GE;
185
186	res.storage.nnz = mat.nonZeros();
187	res.storage.values = mat.derived().valuePtr();
188	res.storage.innerInd = mat.derived().innerIndexPtr();
189	res.storage.outerInd = mat.derived().outerIndexPtr();
190
191	res.setScalarType<typename MatrixType::Scalar>();
192
193	// FIXME the following is not very accurate
194	if (MatrixType::Flags & Upper)
195	res.Mtype = SLU_TRU;
196	if (MatrixType::Flags & Lower)
197	res.Mtype = SLU_TRL;
198
199	eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
200
201	return res;
202	}
203	};
204
205	template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
206	struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
207	{
208	typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
209	static void run(MatrixType& mat, SluMatrix& res)
210	{
211	eigen_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
212	res.setStorageType(SLU_DN);
213	res.setScalarType<Scalar>();
214	res.Mtype = SLU_GE;
215
216	res.nrow = mat.rows();
217	res.ncol = mat.cols();
218
219	res.storage.lda = mat.outerStride();
220	res.storage.values = mat.data();
221	}
222	};
223
224	template<typename Derived>
225	struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
226	{
227	typedef Derived MatrixType;
228	static void run(MatrixType& mat, SluMatrix& res)
229	{
230	if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
231	{
232	res.setStorageType(SLU_NR);
233	res.nrow = mat.cols();
234	res.ncol = mat.rows();
235	}
236	else
237	{
238	res.setStorageType(SLU_NC);
239	res.nrow = mat.rows();
240	res.ncol = mat.cols();
241	}
242
243	res.Mtype = SLU_GE;
244
245	res.storage.nnz = mat.nonZeros();
246	res.storage.values = mat.valuePtr();
247	res.storage.innerInd = mat.innerIndexPtr();
248	res.storage.outerInd = mat.outerIndexPtr();
249
250	res.setScalarType<typename MatrixType::Scalar>();
251
252	// FIXME the following is not very accurate
253	if (MatrixType::Flags & Upper)
254	res.Mtype = SLU_TRU;
255	if (MatrixType::Flags & Lower)
256	res.Mtype = SLU_TRL;
257
258	eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
259	}
260	};
261
262	namespace internal {
263
264	template<typename MatrixType>
265	SluMatrix asSluMatrix(MatrixType& mat)
266	{
267	return SluMatrix::Map(mat);
268	}
269
270	/** View a Super LU matrix as an Eigen expression */
271	template<typename Scalar, int Flags, typename Index>
272	MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat)
273	{
274	eigen_assert((Flags&RowMajor)==RowMajor && sluMat.Stype == SLU_NR
275	\|\| (Flags&ColMajor)==ColMajor && sluMat.Stype == SLU_NC);
276
277	Index outerSize = (Flags&RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow;
278
279	return MappedSparseMatrix<Scalar,Flags,Index>(
280	sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize],
281	sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast<Scalar*>(sluMat.storage.values) );
282	}
283
284	} // end namespace internal
285
286	/** \ingroup SuperLUSupport_Module
287	* \class SuperLUBase
288	* \brief The base class for the direct and incomplete LU factorization of SuperLU
289	*/
290	template<typename _MatrixType, typename Derived>
291	class SuperLUBase : internal::noncopyable
292	{
293	public:
294	typedef _MatrixType MatrixType;
295	typedef typename MatrixType::Scalar Scalar;
296	typedef typename MatrixType::RealScalar RealScalar;
297	typedef typename MatrixType::Index Index;
298	typedef Matrix<Scalar,Dynamic,1> Vector;
299	typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
300	typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
301	typedef SparseMatrix<Scalar> LUMatrixType;
302
303	public:
304
305	SuperLUBase() {}
306
307	~SuperLUBase()
308	{
309	clearFactors();
310	}
311
312	Derived& derived() { return static_cast<Derived>(this); }
313	const Derived& derived() const { return static_cast<const Derived>(this); }
314
315	inline Index rows() const { return m_matrix.rows(); }
316	inline Index cols() const { return m_matrix.cols(); }
317
318	/** \returns a reference to the Super LU option object to configure the Super LU algorithms. */
319	inline superlu_options_t& options() { return m_sluOptions; }
320
321	/** \brief Reports whether previous computation was successful.
322	*
323	* \returns \c Success if computation was succesful,
324	* \c NumericalIssue if the matrix.appears to be negative.
325	*/
326	ComputationInfo info() const
327	{
328	eigen_assert(m_isInitialized && "Decomposition is not initialized.");
329	return m_info;
330	}
331
332	/** Computes the sparse Cholesky decomposition of \a matrix */
333	void compute(const MatrixType& matrix)
334	{
335	derived().analyzePattern(matrix);
336	derived().factorize(matrix);
337	}
338
339	/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
340	*
341	* \sa compute()
342	*/
343	template<typename Rhs>
344	inline const internal::solve_retval<SuperLUBase, Rhs> solve(const MatrixBase<Rhs>& b) const
345	{
346	eigen_assert(m_isInitialized && "SuperLU is not initialized.");
347	eigen_assert(rows()==b.rows()
348	&& "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
349	return internal::solve_retval<SuperLUBase, Rhs>(*this, b.derived());
350	}
351
352	/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
353	*
354	* \sa compute()
355	*/
356	// template<typename Rhs>
357	// inline const internal::sparse_solve_retval<SuperLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const
358	// {
359	// eigen_assert(m_isInitialized && "SuperLU is not initialized.");
360	// eigen_assert(rows()==b.rows()
361	// && "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
362	// return internal::sparse_solve_retval<SuperLU, Rhs>(*this, b.derived());
363	// }
364
365	/** Performs a symbolic decomposition on the sparcity of \a matrix.
366	*
367	* This function is particularly useful when solving for several problems having the same structure.
368	*
369	* \sa factorize()
370	*/
371	void analyzePattern(const MatrixType& /matrix/)
372	{
373	m_isInitialized = true;
374	m_info = Success;
375	m_analysisIsOk = true;
376	m_factorizationIsOk = false;
377	}
378
379	template<typename Stream>
380	void dumpMemory(Stream& s)
381	{}
382
383	protected:
384
385	void initFactorization(const MatrixType& a)
386	{
387	set_default_options(&this->m_sluOptions);
388
389	const int size = a.rows();
390	m_matrix = a;
391
392	m_sluA = internal::asSluMatrix(m_matrix);
393	clearFactors();
394
395	m_p.resize(size);
396	m_q.resize(size);
397	m_sluRscale.resize(size);
398	m_sluCscale.resize(size);
399	m_sluEtree.resize(size);
400
401	// set empty B and X
402	m_sluB.setStorageType(SLU_DN);
403	m_sluB.setScalarType<Scalar>();
404	m_sluB.Mtype = SLU_GE;
405	m_sluB.storage.values = 0;
406	m_sluB.nrow = 0;
407	m_sluB.ncol = 0;
408	m_sluB.storage.lda = size;
409	m_sluX = m_sluB;
410
411	m_extractedDataAreDirty = true;
412	}
413
414	void init()
415	{
416	m_info = InvalidInput;
417	m_isInitialized = false;
418	m_sluL.Store = 0;
419	m_sluU.Store = 0;
420	}
421
422	void extractData() const;
423
424	void clearFactors()
425	{
426	if(m_sluL.Store)
427	Destroy_SuperNode_Matrix(&m_sluL);
428	if(m_sluU.Store)
429	Destroy_CompCol_Matrix(&m_sluU);
430
431	m_sluL.Store = 0;
432	m_sluU.Store = 0;
433
434	memset(&m_sluL,0,sizeof m_sluL);
435	memset(&m_sluU,0,sizeof m_sluU);
436	}
437
438	// cached data to reduce reallocation, etc.
439	mutable LUMatrixType m_l;
440	mutable LUMatrixType m_u;
441	mutable IntColVectorType m_p;
442	mutable IntRowVectorType m_q;
443
444	mutable LUMatrixType m_matrix; // copy of the factorized matrix
445	mutable SluMatrix m_sluA;
446	mutable SuperMatrix m_sluL, m_sluU;
447	mutable SluMatrix m_sluB, m_sluX;
448	mutable SuperLUStat_t m_sluStat;
449	mutable superlu_options_t m_sluOptions;
450	mutable std::vector<int> m_sluEtree;
451	mutable Matrix<RealScalar,Dynamic,1> m_sluRscale, m_sluCscale;
452	mutable Matrix<RealScalar,Dynamic,1> m_sluFerr, m_sluBerr;
453	mutable char m_sluEqued;
454
455	mutable ComputationInfo m_info;
456	bool m_isInitialized;
457	int m_factorizationIsOk;
458	int m_analysisIsOk;
459	mutable bool m_extractedDataAreDirty;
460
461	private:
462	SuperLUBase(SuperLUBase& ) { }
463	};
464
465
466	/** \ingroup SuperLUSupport_Module
467	* \class SuperLU
468	* \brief A sparse direct LU factorization and solver based on the SuperLU library
469	*
470	* This class allows to solve for A.X = B sparse linear problems via a direct LU factorization
471	* using the SuperLU library. The sparse matrix A must be squared and invertible. The vectors or matrices
472	* X and B can be either dense or sparse.
473	*
474	* \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
475	*
476	* \sa \ref TutorialSparseDirectSolvers
477	*/
478	template<typename _MatrixType>
479	class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> >
480	{
481	public:
482	typedef SuperLUBase<_MatrixType,SuperLU> Base;
483	typedef _MatrixType MatrixType;
484	typedef typename Base::Scalar Scalar;
485	typedef typename Base::RealScalar RealScalar;
486	typedef typename Base::Index Index;
487	typedef typename Base::IntRowVectorType IntRowVectorType;
488	typedef typename Base::IntColVectorType IntColVectorType;
489	typedef typename Base::LUMatrixType LUMatrixType;
490	typedef TriangularView<LUMatrixType, Lower\|UnitDiag> LMatrixType;
491	typedef TriangularView<LUMatrixType, Upper> UMatrixType;
492
493	public:
494
495	SuperLU() : Base() { init(); }
496
497	SuperLU(const MatrixType& matrix) : Base()
498	{
499	init();
500	Base::compute(matrix);
501	}
502
503	~SuperLU()
504	{
505	}
506
507	/** Performs a symbolic decomposition on the sparcity of \a matrix.
508	*
509	* This function is particularly useful when solving for several problems having the same structure.
510	*
511	* \sa factorize()
512	*/
513	void analyzePattern(const MatrixType& matrix)
514	{
515	m_info = InvalidInput;
516	m_isInitialized = false;
517	Base::analyzePattern(matrix);
518	}
519
520	/** Performs a numeric decomposition of \a matrix
521	*
522	* The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
523	*
524	* \sa analyzePattern()
525	*/
526	void factorize(const MatrixType& matrix);
527
528	#ifndef EIGEN_PARSED_BY_DOXYGEN
529	/** \internal */
530	template<typename Rhs,typename Dest>
531	void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
532	#endif // EIGEN_PARSED_BY_DOXYGEN
533
534	inline const LMatrixType& matrixL() const
535	{
536	if (m_extractedDataAreDirty) this->extractData();
537	return m_l;
538	}
539
540	inline const UMatrixType& matrixU() const
541	{
542	if (m_extractedDataAreDirty) this->extractData();
543	return m_u;
544	}
545
546	inline const IntColVectorType& permutationP() const
547	{
548	if (m_extractedDataAreDirty) this->extractData();
549	return m_p;
550	}
551
552	inline const IntRowVectorType& permutationQ() const
553	{
554	if (m_extractedDataAreDirty) this->extractData();
555	return m_q;
556	}
557
558	Scalar determinant() const;
559
560	protected:
561
562	using Base::m_matrix;
563	using Base::m_sluOptions;
564	using Base::m_sluA;
565	using Base::m_sluB;
566	using Base::m_sluX;
567	using Base::m_p;
568	using Base::m_q;
569	using Base::m_sluEtree;
570	using Base::m_sluEqued;
571	using Base::m_sluRscale;
572	using Base::m_sluCscale;
573	using Base::m_sluL;
574	using Base::m_sluU;
575	using Base::m_sluStat;
576	using Base::m_sluFerr;
577	using Base::m_sluBerr;
578	using Base::m_l;
579	using Base::m_u;
580
581	using Base::m_analysisIsOk;
582	using Base::m_factorizationIsOk;
583	using Base::m_extractedDataAreDirty;
584	using Base::m_isInitialized;
585	using Base::m_info;
586
587	void init()
588	{
589	Base::init();
590
591	set_default_options(&this->m_sluOptions);
592	m_sluOptions.PrintStat = NO;
593	m_sluOptions.ConditionNumber = NO;
594	m_sluOptions.Trans = NOTRANS;
595	m_sluOptions.ColPerm = COLAMD;
596	}
597
598
599	private:
600	SuperLU(SuperLU& ) { }
601	};
602
603	template<typename MatrixType>
604	void SuperLU<MatrixType>::factorize(const MatrixType& a)
605	{
606	eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
607	if(!m_analysisIsOk)
608	{
609	m_info = InvalidInput;
610	return;
611	}
612
613	this->initFactorization(a);
614
615	int info = 0;
616	RealScalar recip_pivot_growth, rcond;
617	RealScalar ferr, berr;
618
619	StatInit(&m_sluStat);
620	SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
621	&m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
622	&m_sluL, &m_sluU,
623	NULL, 0,
624	&m_sluB, &m_sluX,
625	&recip_pivot_growth, &rcond,
626	&ferr, &berr,
627	&m_sluStat, &info, Scalar());
628	StatFree(&m_sluStat);
629
630	m_extractedDataAreDirty = true;
631
632	// FIXME how to better check for errors ???
633	m_info = info == 0 ? Success : NumericalIssue;
634	m_factorizationIsOk = true;
635	}
636
637	template<typename MatrixType>
638	template<typename Rhs,typename Dest>
639	void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
640	{
641	eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
642
643	const int size = m_matrix.rows();
644	const int rhsCols = b.cols();
645	eigen_assert(size==b.rows());
646
647	m_sluOptions.Trans = NOTRANS;
648	m_sluOptions.Fact = FACTORED;
649	m_sluOptions.IterRefine = NOREFINE;
650
651
652	m_sluFerr.resize(rhsCols);
653	m_sluBerr.resize(rhsCols);
654	m_sluB = SluMatrix::Map(b.const_cast_derived());
655	m_sluX = SluMatrix::Map(x.derived());
656
657	typename Rhs::PlainObject b_cpy;
658	if(m_sluEqued!='N')
659	{
660	b_cpy = b;
661	m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
662	}
663
664	StatInit(&m_sluStat);
665	int info = 0;
666	RealScalar recip_pivot_growth, rcond;
667	SuperLU_gssvx(&m_sluOptions, &m_sluA,
668	m_q.data(), m_p.data(),
669	&m_sluEtree[0], &m_sluEqued,
670	&m_sluRscale[0], &m_sluCscale[0],
671	&m_sluL, &m_sluU,
672	NULL, 0,
673	&m_sluB, &m_sluX,
674	&recip_pivot_growth, &rcond,
675	&m_sluFerr[0], &m_sluBerr[0],
676	&m_sluStat, &info, Scalar());
677	StatFree(&m_sluStat);
678	m_info = info==0 ? Success : NumericalIssue;
679	}
680
681	// the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
682	//
683	// Copyright (c) 1994 by Xerox Corporation. All rights reserved.
684	//
685	// THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
686	// EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
687	//
688	template<typename MatrixType, typename Derived>
689	void SuperLUBase<MatrixType,Derived>::extractData() const
690	{
691	eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
692	if (m_extractedDataAreDirty)
693	{
694	int upper;
695	int fsupc, istart, nsupr;
696	int lastl = 0, lastu = 0;
697	SCformat Lstore = static_cast<SCformat>(m_sluL.Store);
698	NCformat Ustore = static_cast<NCformat>(m_sluU.Store);
699	Scalar *SNptr;
700
701	const int size = m_matrix.rows();
702	m_l.resize(size,size);
703	m_l.resizeNonZeros(Lstore->nnz);
704	m_u.resize(size,size);
705	m_u.resizeNonZeros(Ustore->nnz);
706
707	int* Lcol = m_l.outerIndexPtr();
708	int* Lrow = m_l.innerIndexPtr();
709	Scalar* Lval = m_l.valuePtr();
710
711	int* Ucol = m_u.outerIndexPtr();
712	int* Urow = m_u.innerIndexPtr();
713	Scalar* Uval = m_u.valuePtr();
714
715	Ucol[0] = 0;
716	Ucol[0] = 0;
717
718	/* for each supernode */
719	for (int k = 0; k <= Lstore->nsuper; ++k)
720	{
721	fsupc = L_FST_SUPC(k);
722	istart = L_SUB_START(fsupc);
723	nsupr = L_SUB_START(fsupc+1) - istart;
724	upper = 1;
725
726	/* for each column in the supernode */
727	for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
728	{
729	SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
730
731	/* Extract U */
732	for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
733	{
734	Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
735	/* Matlab doesn't like explicit zero. */
736	if (Uval[lastu] != 0.0)
737	Urow[lastu++] = U_SUB(i);
738	}
739	for (int i = 0; i < upper; ++i)
740	{
741	/* upper triangle in the supernode */
742	Uval[lastu] = SNptr[i];
743	/* Matlab doesn't like explicit zero. */
744	if (Uval[lastu] != 0.0)
745	Urow[lastu++] = L_SUB(istart+i);
746	}
747	Ucol[j+1] = lastu;
748
749	/* Extract L */
750	Lval[lastl] = 1.0; /* unit diagonal */
751	Lrow[lastl++] = L_SUB(istart + upper - 1);
752	for (int i = upper; i < nsupr; ++i)
753	{
754	Lval[lastl] = SNptr[i];
755	/* Matlab doesn't like explicit zero. */
756	if (Lval[lastl] != 0.0)
757	Lrow[lastl++] = L_SUB(istart+i);
758	}
759	Lcol[j+1] = lastl;
760
761	++upper;
762	} /* for j ... */
763
764	} /* for k ... */
765
766	// squeeze the matrices :
767	m_l.resizeNonZeros(lastl);
768	m_u.resizeNonZeros(lastu);
769
770	m_extractedDataAreDirty = false;
771	}
772	}
773
774	template<typename MatrixType>
775	typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const
776	{
777	eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for computing the determinant, you must first call either compute() or analyzePattern()/factorize()");
778
779	if (m_extractedDataAreDirty)
780	this->extractData();
781
782	Scalar det = Scalar(1);
783	for (int j=0; j<m_u.cols(); ++j)
784	{
785	if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0)
786	{
787	int lastId = m_u.outerIndexPtr()[j+1]-1;
788	eigen_assert(m_u.innerIndexPtr()[lastId]<=j);
789	if (m_u.innerIndexPtr()[lastId]==j)
790	det *= m_u.valuePtr()[lastId];
791	}
792	}
793	if(m_sluEqued!='N')
794	return det/m_sluRscale.prod()/m_sluCscale.prod();
795	else
796	return det;
797	}
798
799	#ifdef EIGEN_PARSED_BY_DOXYGEN
800	#define EIGEN_SUPERLU_HAS_ILU
801	#endif
802
803	#ifdef EIGEN_SUPERLU_HAS_ILU
804
805	/** \ingroup SuperLUSupport_Module
806	* \class SuperILU
807	* \brief A sparse direct \b incomplete LU factorization and solver based on the SuperLU library
808	*
809	* This class allows to solve for an approximate solution of A.X = B sparse linear problems via an incomplete LU factorization
810	* using the SuperLU library. This class is aimed to be used as a preconditioner of the iterative linear solvers.
811	*
812	* \warning This class requires SuperLU 4 or later.
813	*
814	* \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
815	*
816	* \sa \ref TutorialSparseDirectSolvers, class ConjugateGradient, class BiCGSTAB
817	*/
818
819	template<typename _MatrixType>
820	class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
821	{
822	public:
823	typedef SuperLUBase<_MatrixType,SuperILU> Base;
824	typedef _MatrixType MatrixType;
825	typedef typename Base::Scalar Scalar;
826	typedef typename Base::RealScalar RealScalar;
827	typedef typename Base::Index Index;
828
829	public:
830
831	SuperILU() : Base() { init(); }
832
833	SuperILU(const MatrixType& matrix) : Base()
834	{
835	init();
836	Base::compute(matrix);
837	}
838
839	~SuperILU()
840	{
841	}
842
843	/** Performs a symbolic decomposition on the sparcity of \a matrix.
844	*
845	* This function is particularly useful when solving for several problems having the same structure.
846	*
847	* \sa factorize()
848	*/
849	void analyzePattern(const MatrixType& matrix)
850	{
851	Base::analyzePattern(matrix);
852	}
853
854	/** Performs a numeric decomposition of \a matrix
855	*
856	* The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
857	*
858	* \sa analyzePattern()
859	*/
860	void factorize(const MatrixType& matrix);
861
862	#ifndef EIGEN_PARSED_BY_DOXYGEN
863	/** \internal */
864	template<typename Rhs,typename Dest>
865	void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
866	#endif // EIGEN_PARSED_BY_DOXYGEN
867
868	protected:
869
870	using Base::m_matrix;
871	using Base::m_sluOptions;
872	using Base::m_sluA;
873	using Base::m_sluB;
874	using Base::m_sluX;
875	using Base::m_p;
876	using Base::m_q;
877	using Base::m_sluEtree;
878	using Base::m_sluEqued;
879	using Base::m_sluRscale;
880	using Base::m_sluCscale;
881	using Base::m_sluL;
882	using Base::m_sluU;
883	using Base::m_sluStat;
884	using Base::m_sluFerr;
885	using Base::m_sluBerr;
886	using Base::m_l;
887	using Base::m_u;
888
889	using Base::m_analysisIsOk;
890	using Base::m_factorizationIsOk;
891	using Base::m_extractedDataAreDirty;
892	using Base::m_isInitialized;
893	using Base::m_info;
894
895	void init()
896	{
897	Base::init();
898
899	ilu_set_default_options(&m_sluOptions);
900	m_sluOptions.PrintStat = NO;
901	m_sluOptions.ConditionNumber = NO;
902	m_sluOptions.Trans = NOTRANS;
903	m_sluOptions.ColPerm = MMD_AT_PLUS_A;
904
905	// no attempt to preserve column sum
906	m_sluOptions.ILU_MILU = SILU;
907	// only basic ILU(k) support -- no direct control over memory consumption
908	// better to use ILU_DropRule = DROP_BASIC \| DROP_AREA
909	// and set ILU_FillFactor to max memory growth
910	m_sluOptions.ILU_DropRule = DROP_BASIC;
911	m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10;
912	}
913
914	private:
915	SuperILU(SuperILU& ) { }
916	};
917
918	template<typename MatrixType>
919	void SuperILU<MatrixType>::factorize(const MatrixType& a)
920	{
921	eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
922	if(!m_analysisIsOk)
923	{
924	m_info = InvalidInput;
925	return;
926	}
927
928	this->initFactorization(a);
929
930	int info = 0;
931	RealScalar recip_pivot_growth, rcond;
932
933	StatInit(&m_sluStat);
934	SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
935	&m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
936	&m_sluL, &m_sluU,
937	NULL, 0,
938	&m_sluB, &m_sluX,
939	&recip_pivot_growth, &rcond,
940	&m_sluStat, &info, Scalar());
941	StatFree(&m_sluStat);
942
943	// FIXME how to better check for errors ???
944	m_info = info == 0 ? Success : NumericalIssue;
945	m_factorizationIsOk = true;
946	}
947
948	template<typename MatrixType>
949	template<typename Rhs,typename Dest>
950	void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
951	{
952	eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
953
954	const int size = m_matrix.rows();
955	const int rhsCols = b.cols();
956	eigen_assert(size==b.rows());
957
958	m_sluOptions.Trans = NOTRANS;
959	m_sluOptions.Fact = FACTORED;
960	m_sluOptions.IterRefine = NOREFINE;
961
962	m_sluFerr.resize(rhsCols);
963	m_sluBerr.resize(rhsCols);
964	m_sluB = SluMatrix::Map(b.const_cast_derived());
965	m_sluX = SluMatrix::Map(x.derived());
966
967	typename Rhs::PlainObject b_cpy;
968	if(m_sluEqued!='N')
969	{
970	b_cpy = b;
971	m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
972	}
973
974	int info = 0;
975	RealScalar recip_pivot_growth, rcond;
976
977	StatInit(&m_sluStat);
978	SuperLU_gsisx(&m_sluOptions, &m_sluA,
979	m_q.data(), m_p.data(),
980	&m_sluEtree[0], &m_sluEqued,
981	&m_sluRscale[0], &m_sluCscale[0],
982	&m_sluL, &m_sluU,
983	NULL, 0,
984	&m_sluB, &m_sluX,
985	&recip_pivot_growth, &rcond,
986	&m_sluStat, &info, Scalar());
987	StatFree(&m_sluStat);
988
989	m_info = info==0 ? Success : NumericalIssue;
990	}
991	#endif
992
993	namespace internal {
994
995	template<typename _MatrixType, typename Derived, typename Rhs>
996	struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
997	: solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
998	{
999	typedef SuperLUBase<_MatrixType,Derived> Dec;
1000	EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
1001
1002	template<typename Dest> void evalTo(Dest& dst) const
1003	{
1004	dec().derived()._solve(rhs(),dst);
1005	}
1006	};
1007
1008	template<typename _MatrixType, typename Derived, typename Rhs>
1009	struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
1010	: sparse_solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
1011	{
1012	typedef SuperLUBase<_MatrixType,Derived> Dec;
1013	EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
1014
1015	template<typename Dest> void evalTo(Dest& dst) const
1016	{
1017	dec().derived()._solve(rhs(),dst);
1018	}
1019	};
1020
1021	} // end namespace internal
1022
1023	} // end namespace Eigen
1024
1025	#endif // EIGEN_SUPERLUSUPPORT_H

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Update cookies preferences