Rheolef  7.1
an efficient C++ finite element environment
solver_trilinos_ifpack.cc
Go to the documentation of this file.
1 // 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)
26 #include "solver_trilinos_ifpack.h"
27 #include "Ifpack.h"
28 #include "Epetra_Vector.h"
29 
30 namespace rheolef {
31 using namespace std;
32 
33 // convert from rheolef to trilinos
34 static
35 Epetra_CrsMatrix*
36 csr2petra (const csr<double,sequential>& a, Epetra_Map*& petra_ownership_ptr)
37 {
38  error_macro ("not yet");
39  return 0;
40 }
41 static
42 Epetra_CrsMatrix*
43 csr2petra (const csr<double,distributed>& a, Epetra_Map*& petra_ownership_ptr)
44 {
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 }
83 template<class T, class M>
84 void
85 solver_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 }
90 template<class T, class M>
91 void
92 solver_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 }
130 template<class T, class M>
131 solver_trilinos_ifpack_rep<T,M>::~solver_trilinos_ifpack_rep ()
132 {
133  destroy_values ();
134 }
135 template<class T, class M>
136 solver_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 }
144 template<class T, class M>
145 void
146 solver_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 }
157 template<class T, class M>
158 vec<T,M>
159 solver_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 }
165 template<class T, class M>
166 vec<T,M>
167 solver_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 
178 template class solver_trilinos_ifpack_rep<double,sequential>;
179 #ifdef _RHEOLEF_HAVE_MPI
180 template 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:425
rheolef::std type
size_t size_type
Definition: basis_get.cc:76
Expr1::float_type T
Definition: field_expr.h:261
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
Definition: sphere.icc:25
Expr1::memory_type M
Definition: vec_expr_v2.h:416