24#include "rheolef/config.h"
25#if defined(_RHEOLEF_HAVE_TRILINOS) && defined(_RHEOLEF_HAVE_MPI)
28#include "Epetra_Vector.h"
36csr2petra (
const csr<double,sequential>& a, Epetra_Map*& petra_ownership_ptr)
43csr2petra (
const csr<double,distributed>& a, Epetra_Map*& petra_ownership_ptr)
45 typedef csr<double,distributed>::size_type
size_type;
46 distributor ownership =
a.row_ownership();
47 Epetra_MpiComm petra_comm (ownership.comm());
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;
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]);
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();
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;
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);
77 a_petra.InsertGlobalValues(dis_i, nz_by_row[i], val.begin().operator->(), jdx.begin().operator->());
80 a_petra.FillComplete(optimize);
83template<
class T,
class M>
85solver_trilinos_ifpack_rep<T,M>::destroy_values ()
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; }
90template<
class T,
class M>
92solver_trilinos_ifpack_rep<T,M>::update_values (
const csr<T,M>& a)
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);
108 string type = (
a.is_symmetric() ?
"IC" :
"ILU");
110 Epetra_CrsMatrix* a_petra_ptr = csr2petra (a, _petra_ownership_ptr);
111 _ilu_ptr = factory.Create (type, a_petra_ptr);
112 Teuchos::ParameterList params;
122 params.set(
"fact: level-of-fill", 0);
123 _ilu_ptr -> SetParameters (params);
124 _ilu_ptr -> Initialize();
125 _ilu_ptr -> Compute();
127 delete_macro(a_petra_ptr);
130template<
class T,
class M>
131solver_trilinos_ifpack_rep<T,M>::~solver_trilinos_ifpack_rep ()
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),
139 _petra_ownership_ptr(0)
144template<
class T,
class M>
146solver_trilinos_ifpack_rep<T,M>::solve (
const vec<T,M>& b,
bool do_transpose, vec<T,M>& x)
const
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);
157template<
class T,
class M>
159solver_trilinos_ifpack_rep<T,M>::solve (
const vec<T,M>& b)
const
161 vec<T,M> x (
b.ownership());
165template<
class T,
class M>
167solver_trilinos_ifpack_rep<T,M>::trans_solve (
const vec<T,M>& b)
const
169 vec<T,M> x (
b.ownership());
178template class solver_trilinos_ifpack_rep<double,sequential>;
179#ifdef _RHEOLEF_HAVE_MPI
180template class solver_trilinos_ifpack_rep<double,distributed>;
field::size_type size_type
#define error_macro(message)
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)