Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
solver_trilinos_ifpack.cc
Go to the documentation of this file.
1
21// preconditioner based on trilinos/ifpack, both seq & mpi implementations
22// note: the library requires mpi to be present => no seq-only implementation
23//
24#include "rheolef/config.h"
25#if defined(_RHEOLEF_HAVE_TRILINOS) && defined(_RHEOLEF_HAVE_MPI)
27#include "Ifpack.h"
28#include "Epetra_Vector.h"
29
30namespace rheolef {
31using namespace std;
32
33// convert from rheolef to trilinos
34static
35Epetra_CrsMatrix*
36csr2petra (const csr<double,sequential>& a, Epetra_Map*& petra_ownership_ptr)
37{
38 error_macro ("not yet");
39 return 0;
40}
41static
42Epetra_CrsMatrix*
43csr2petra (const csr<double,distributed>& a, Epetra_Map*& petra_ownership_ptr)
44{
45 typedef csr<double,distributed>::size_type size_type;
46 distributor ownership = a.row_ownership();
47 Epetra_MpiComm petra_comm (ownership.comm()); // or Epetra_SerialComm
48 petra_ownership_ptr = new_macro (Epetra_Map (ownership.dis_size(), ownership.size(), 0, petra_comm));
49 std::vector<int> nz_by_row (a.nrow());
50 csr<double,distributed>::const_iterator dia_ia = a.begin();
51 csr<double,distributed>::const_iterator ext_ia = a.ext_begin();
52 int max_nz_by_row = 0;
53 int nnz = 0;
54 for (size_type i = 0, n = a.nrow(); i < n; i++) {
55 nz_by_row[i] = (dia_ia[i+1] - dia_ia[i]) + (ext_ia[i+1] - ext_ia[i]);
56 max_nz_by_row = std::max(max_nz_by_row, nz_by_row[i]);
57 nnz += nz_by_row[i];
58 }
59 bool static_storage = true;
60 Epetra_CrsMatrix* a_petra_ptr = new_macro (Epetra_CrsMatrix (Copy, *petra_ownership_ptr, nz_by_row.begin().operator->(), static_storage));
61 Epetra_CrsMatrix& a_petra = *a_petra_ptr;
62 std::vector<int> jdx (max_nz_by_row);
63 std::vector<double> val (max_nz_by_row);
64 size_type first_i = ownership.first_index();
65 for (size_type i = 0, n = a.nrow(); i < n; i++) {
66 size_type q = 0;
67 for (csr<double,distributed>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++, q++) {
68 jdx[q] = (*p).first + first_i;
69 val[q] = (*p).second;
70 }
71 for (csr<double,distributed>::const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++, q++) {
72 jdx[q] = a.jext2dis_j ((*p).first);
73 val[q] = (*p).second;
74 }
75 size_type dis_i = first_i + i;
76 check_macro (int(q) == nz_by_row[i], "unexpected");
77 a_petra.InsertGlobalValues(dis_i, nz_by_row[i], val.begin().operator->(), jdx.begin().operator->());
78 }
79 bool optimize = true;
80 a_petra.FillComplete(optimize);
81 return a_petra_ptr;
82}
83template<class T, class M>
84void
85solver_trilinos_ifpack_rep<T,M>::destroy_values ()
86{
87 if (_ilu_ptr) { delete_macro(_ilu_ptr); _ilu_ptr = 0; }
88 if (_petra_ownership_ptr) { delete_macro(_petra_ownership_ptr); _petra_ownership_ptr = 0; }
89}
90template<class T, class M>
91void
92solver_trilinos_ifpack_rep<T,M>::update_values (const csr<T,M>& a)
93{
94 destroy_values ();
95#ifdef TO_CLEAN
96 // ILUT, with thresold, cannot exploit the constant sparsity pattern of a
97 // thus, the factorization may be recomputed completely
98 double k = 1; // TODO: in options
99 double drop_tol = 0;
100 int nnz = a.dis_nnz();
101 int n = a.dis_nrow();
102 int level_fill = int (k*nnz/(2.*n) + 1);
103 check_macro (a.is_symmetric(), "ict(k,e): unsymmetric matrix not supported");
104 _ilu_ptr = new_macro (Ifpack_CrsIct (*a_petra_ptr, drop_tol, level_fill));
105 _ilu_ptr->InitValues(*a_petra_ptr);
106 _ilu_ptr->Factor();
107#endif // TO_CLEAN
108 string type = (a.is_symmetric() ? "IC" : "ILU"); // TODO: ICT & ILUT
109 Ifpack factory;
110 Epetra_CrsMatrix* a_petra_ptr = csr2petra (a, _petra_ownership_ptr);
111 _ilu_ptr = factory.Create (type, a_petra_ptr);
112 Teuchos::ParameterList params;
113 // TODO: fact parameters
114 // TODO: overlap with domains
115 // fact: level-of-fill for IC and ILU.
116 // fact: ict level-of-fill [double] Level-of-fill for ICT.
117 // fact: ilut level-of-fill [double] Level-of-fill for ILUT.
118 // fact: relax value [double] Relaxation value.
119 // fact: absolute threshold [double] Value of alpha in equation (16).
120 // fact: level-of-fill relative threshold [double] Value of rho in equation (16).
121 // B = alpha*sgn(A) + rho*A (16)
122 params.set("fact: level-of-fill", 0);
123 _ilu_ptr -> SetParameters (params);
124 _ilu_ptr -> Initialize();
125 _ilu_ptr -> Compute();
126#ifdef TO_CLEAN
127 delete_macro(a_petra_ptr); // deleted by _ilu_ptr ?
128#endif // TO_CLEAN
129}
130template<class T, class M>
131solver_trilinos_ifpack_rep<T,M>::~solver_trilinos_ifpack_rep ()
132{
133 destroy_values ();
134}
135template<class T, class M>
136solver_trilinos_ifpack_rep<T,M>::solver_trilinos_ifpack_rep (const csr<T,M>& a, const solver_option& opt)
137 : solver_abstract_rep<T,M>(opt),
138 _ilu_ptr(0),
139 _petra_ownership_ptr(0)
140{
141 // TODO: (k, drop_tol) in options
142 update_values (a);
143}
144template<class T, class M>
145void
146solver_trilinos_ifpack_rep<T,M>::solve (const vec<T,M>& b, bool do_transpose, vec<T,M>& x) const
147{
148 const double* b_values = b.begin().operator->();
149 double* x_values = x.begin().operator->();
150 check_macro (int(x.ownership().size()) == _petra_ownership_ptr -> NumMyElements(),
151 "incomplete choleski preconditioner: incompatible right-hand size");
152 Epetra_Vector b_petra (View, *_petra_ownership_ptr, const_cast<double*>(b_values));
153 Epetra_Vector x_petra (View, *_petra_ownership_ptr, x_values);
154 _ilu_ptr -> SetUseTranspose (do_transpose);
155 _ilu_ptr -> ApplyInverse (b_petra, x_petra);
156}
157template<class T, class M>
158vec<T,M>
159solver_trilinos_ifpack_rep<T,M>::solve (const vec<T,M>& b) const
160{
161 vec<T,M> x (b.ownership());
162 solve (b, false, x);
163 return x;
164}
165template<class T, class M>
166vec<T,M>
167solver_trilinos_ifpack_rep<T,M>::trans_solve (const vec<T,M>& b) const
168{
169 vec<T,M> x (b.ownership());
170 solve (b, true, x);
171 return x;
172}
173// ----------------------------------------------------------------------------
174// instanciation in library
175// ----------------------------------------------------------------------------
176// TODO: code is only valid here for T=double
177
178template class solver_trilinos_ifpack_rep<double,sequential>;
179#ifdef _RHEOLEF_HAVE_MPI
180template class solver_trilinos_ifpack_rep<double,distributed>;
181#endif // _RHEOLEF_HAVE_MPI
182
183} // namespace rheolef
184#endif // _RHEOLEF_HAVE_TRILINOS && _RHEOLEF_HAVE_MPI
field::size_type size_type
Definition branch.cc:430
rheolef::std type
#define error_macro(message)
Definition dis_macros.h:49
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
void solve(tiny_matrix< T > &a, tiny_vector< size_t > &piv, const tiny_vector< T > &b, tiny_vector< T > &x)
Definition tiny_lu.h:92
STL namespace.
Definition sphere.icc:25
Expr1::memory_type M