IlutPreconditioning.hxx
1 // Copyright (C) 2010 Marc Duruflé
2 //
3 // This file is part of the linear-algebra library Seldon,
4 // http://seldon.sourceforge.net/.
5 //
6 // Seldon is free software; you can redistribute it and/or modify it under the
7 // terms of the GNU Lesser General Public License as published by the Free
8 // Software Foundation; either version 2.1 of the License, or (at your option)
9 // any later version.
10 //
11 // Seldon is distributed in the hope that it will be useful, but WITHOUT ANY
12 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
14 // more details.
15 //
16 // You should have received a copy of the GNU Lesser General Public License
17 // along with Seldon. If not, see http://www.gnu.org/licenses/.
18 
19 
20 #ifndef SELDON_FILE_ILUT_PRECONDITIONING_HXX
21 
22 namespace Seldon
23 {
24 
25  template<class T, class Allocator
26  = typename SeldonDefaultAllocator<ArrayRowSparse, T>::allocator>
27  class IlutPreconditioning : public Preconditioner_Base<T>,
28  public VirtualSparseDirectSolver<T>
29  {
30  protected :
32  int print_level;
34  bool symmetric_algorithm;
36  int type_ilu;
40  int fill_level;
44  int additional_fill;
46  int mbloc;
48  typename ClassComplexType<T>::Treal alpha;
50  typename ClassComplexType<T>::Treal droptol;
52  typename ClassComplexType<T>::Treal permtol;
54  IVect permutation_row, permutation_col;
56  Matrix<T, Symmetric, ArrayRowSymSparse, Allocator> mat_sym;
58  Matrix<T, General, ArrayRowSparse, Allocator> mat_unsym;
59 
60  public :
61 
63  enum {ILUT, ILU_D, ILUT_K, ILU_0, MILU_0, ILU_K};
64 
65  IlutPreconditioning();
66 
67  void Clear();
68 
69  bool UseInteger8() const;
70  void HideMessages();
71  void ShowMessages();
72 
73  int GetFactorisationType() const;
74  int GetFillLevel() const;
75  int GetAdditionalFillNumber() const;
76  int GetPrintLevel() const;
77  int GetPivotBlockInteger() const;
78  size_t GetMemorySize() const;
79  int GetInfoFactorization() const;
80 
81  void SetFactorisationType(int);
82  void SetFillLevel(int);
83  void SetAdditionalFillNumber(int);
84  void SetPrintLevel(int);
85  void SetPivotBlockInteger(int);
86  void SetSymmetricAlgorithm();
87  void SetUnsymmetricAlgorithm();
88 
89  typename ClassComplexType<T>::Treal GetDroppingThreshold() const;
90  typename ClassComplexType<T>::Treal GetDiagonalCoefficient() const;
91  typename ClassComplexType<T>::Treal GetPivotThreshold() const;
92 
93  void SetDroppingThreshold(typename ClassComplexType<T>::Treal);
94  void SetDiagonalCoefficient(typename ClassComplexType<T>::Treal);
95  void SetPivotThreshold(double);
96 
97  template<class MatrixSparse>
98  void FactorizeSymMatrix(const IVect& perm,
99  MatrixSparse& mat, bool keep_matrix = false);
100 
101  template<class MatrixSparse>
102  void FactorizeUnsymMatrix(const IVect& perm,
103  MatrixSparse& mat, bool keep_matrix = false);
104 
105  template<class T0, class Storage0, class Allocator0>
106  void FactorizeMatrix(const IVect& perm,
107  Matrix<T0, General, Storage0, Allocator0>& mat,
108  bool keep_matrix = false);
109 
110  template<class T0, class Storage0, class Allocator0>
111  void FactorizeMatrix(const IVect& perm,
112  Matrix<T0, Symmetric, Storage0, Allocator0>& mat,
113  bool keep_matrix = false);
114 
115 #ifdef SELDON_WITH_VIRTUAL
116  void Solve(const VirtualMatrix<T>&, const Vector<T>& r, Vector<T>& z);
117  void TransSolve(const VirtualMatrix<T>&, const Vector<T>& r, Vector<T>&);
118 #else
119  template<class Matrix1, class Vector1>
120  void TransSolve(const Matrix1& A, const Vector1& r, Vector1& z);
121 
122  template<class Matrix1, class Vector1>
123  void Solve(const Matrix1& A, const Vector1& r, Vector1& z);
124 #endif
125 
126  template<class Vector1>
127  void TransSolve(Vector1& z);
128 
129  template<class Vector1>
130  void Solve(Vector1& z);
131 
132  template<class Vector1>
133  void Solve(const SeldonTranspose&, Vector1& z);
134 
135  void Solve(const SeldonTranspose&, T* x_ptr, int nrhs);
136 
137  };
138 
139  template<class real, class cplx, class Storage, class Allocator>
140  void qsplit_ilut(Vector<cplx, Storage, Allocator>& a, IVect& ind, int first,
141  int n, int ncut, const real& abs_ncut);
142 
143  template<class cplx, class Allocator1, class Allocator2>
144  void GetIlut(const IlutPreconditioning<cplx, Allocator1>& param,
145  Matrix<cplx, General, ArrayRowSparse, Allocator2>& A,
146  IVect& iperm, IVect& rperm);
147 
148  template<class cplx, class Allocator>
149  void GetIluk(int lfil, Matrix<cplx, General, ArrayRowSparse, Allocator>& A);
150 
151  template<class cplx, class Allocator>
152  void GetIlu0(Matrix<cplx, General, ArrayRowSparse, Allocator>& A);
153 
154  template<class cplx, class Allocator>
155  void GetMilu0(Matrix<cplx, General, ArrayRowSparse, Allocator>& A);
156 
157  template<class cplx, class Allocator1, class Allocator2>
158  void GetIlut(const IlutPreconditioning<cplx, Allocator1>& param,
159  Matrix<cplx, Symmetric, ArrayRowSymSparse, Allocator2>& A);
160 
161  template<class cplx, class Allocator>
162  void GetIluk(int lfil, int print_level,
163  Matrix<cplx, Symmetric, ArrayRowSymSparse, Allocator>& A);
164 
165  template<class cplx, class Allocator>
166  void GetIlu0(Matrix<cplx, Symmetric, ArrayRowSymSparse, Allocator>& A);
167 
168  template<class cplx, class Allocator>
169  void GetMilu0(Matrix<cplx, Symmetric, ArrayRowSymSparse, Allocator>& A);
170 
171  template<class MatrixSparse, class T, class Alloc2>
172  void GetLU(MatrixSparse& A, IlutPreconditioning<T, Alloc2>& mat_lu,
173  IVect& permut, bool keep_matrix, T& x);
174 
175  template<class MatrixSparse, class T, class Alloc2>
176  void GetLU(MatrixSparse& A, IlutPreconditioning<T, Alloc2>& mat_lu,
177  IVect& permut, bool keep_matrix, complex<T>& x);
178 
179  template<class MatrixSparse, class T, class Alloc2>
180  void GetLU(MatrixSparse& A, IlutPreconditioning<complex<T>, Alloc2>& mat_lu,
181  IVect& permut, bool keep_matrix, T& x);
182 
183  template<class T0, class Prop, class Storage, class Allocator,
184  class T, class Alloc2>
185  void GetLU(Matrix<T0, Prop, Storage, Allocator>& A,
186  IlutPreconditioning<T, Alloc2>& mat_lu,
187  IVect& permut, bool keep_matrix = false);
188 
189 }
190 
191 #define SELDON_FILE_ILUT_PRECONDITIONING_HXX
192 #endif
Seldon::SeldonTranspose
Definition: MatrixFlag.hxx:32
Seldon::Vector< int, VectFull >
Seldon::IlutPreconditioning::Solve
void Solve(const Matrix1 &A, const Vector1 &r, Vector1 &z)
Applies ilut preconditioning.
Definition: IlutPreconditioning.cxx:418
Seldon::Matrix
Definition: SeldonHeader.hxx:226
Seldon::IlutPreconditioning::mat_unsym
Matrix< T, General, ArrayRowSparse, Allocator > mat_unsym
Unsymmetric matrix.
Definition: IlutPreconditioning.hxx:58
Seldon::IlutPreconditioning::permutation_row
IVect permutation_row
Permutation arrays.
Definition: IlutPreconditioning.hxx:54
Seldon::IlutPreconditioning::alpha
ClassComplexType< T >::Treal alpha
Diagonal compensation parameter (alpha = 0 -> ILU, alpha = 1 -> MILU).
Definition: IlutPreconditioning.hxx:48
Seldon::IlutPreconditioning::SetPivotThreshold
void SetPivotThreshold(double)
Sets the threshold for pivot.
Definition: IlutPreconditioning.cxx:203
Seldon::IlutPreconditioning::print_level
int print_level
Verbosity level.
Definition: IlutPreconditioning.hxx:32
Seldon::IlutPreconditioning::permtol
ClassComplexType< T >::Treal permtol
Threshold for pivoting.
Definition: IlutPreconditioning.hxx:52
Seldon::GetIlut
void GetIlut(const IlutPreconditioning< cplx, Allocator1 > &param, Matrix< cplx, General, ArrayRowSparse, Allocator2 > &A, IVect &iperm, IVect &rperm)
Incomplete factorization with pivot for unsymmetric matrix.
Definition: IlutPreconditioning.cxx:611
Seldon::IlutPreconditioning::mat_sym
Matrix< T, Symmetric, ArrayRowSymSparse, Allocator > mat_sym
Symmetric matrix.
Definition: IlutPreconditioning.hxx:56
Seldon::IlutPreconditioning::fill_level
int fill_level
Maximum number of elements on a row of L or U. For Ilu(k), fill_level = k.
Definition: IlutPreconditioning.hxx:40
Seldon::IlutPreconditioning::additional_fill
int additional_fill
Additional number of elements for each row. This number is only used for ILUT(k)
Definition: IlutPreconditioning.hxx:44
Seldon::GetLU
void GetLU(Matrix< T0, Prop0, Storage0, Allocator0 > &A)
Returns the LU factorization of a matrix.
Definition: Functions_Matrix.cxx:2073
Seldon::VirtualSparseDirectSolver
Base class for an interface with a direct solver.
Definition: SparseSolver.hxx:30
Seldon::IlutPreconditioning::droptol
ClassComplexType< T >::Treal droptol
Threshold used for dropping small terms.
Definition: IlutPreconditioning.hxx:50
Seldon::IlutPreconditioning::symmetric_algorithm
bool symmetric_algorithm
True if symmetric matrix is constructed.
Definition: IlutPreconditioning.hxx:34
Seldon::IlutPreconditioning::type_ilu
int type_ilu
Type of incomplete factorization.
Definition: IlutPreconditioning.hxx:36
Seldon::IlutPreconditioning::mbloc
int mbloc
Size of block where the pivot is searched.
Definition: IlutPreconditioning.hxx:46
Seldon::IlutPreconditioning
Definition: IlutPreconditioning.hxx:27
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::IlutPreconditioning::TransSolve
void TransSolve(const Matrix1 &A, const Vector1 &r, Vector1 &z)
Applies transpose of ilut preconditioning.
Definition: IlutPreconditioning.cxx:449
Seldon::VirtualMatrix
Abstract base class for all matrices.
Definition: Matrix_Base.hxx:42
Seldon::Preconditioner_Base
Base class for preconditioners.
Definition: Iterative.hxx:26