Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
solver_pastix.h
Go to the documentation of this file.
1#ifndef _RHEOLEF_SOLVER_PASTIX_H
2#define _RHEOLEF_SOLVER_PASTIX_H
23// solver pastix implementation: interface
24//
25#include "rheolef/config.h"
26
27#ifdef _RHEOLEF_HAVE_PASTIX
28
29#include "rheolef/solver.h"
30extern "C" {
31#define COMPLEXFLOAT_ /* workarroud a compile problem here */
32#define COMPLEXDOUBLE_
33#ifndef _RHEOLEF_HAVE_MPI
34#define FORCE_NOMPI
35#define MPI_Comm int
36#endif // _RHEOLEF_HAVE_MPI
37#include "pastix.h"
38#include "cscd_utils.h"
39#undef COMPLEXFLOAT_
40#undef COMPLEXDOUBLE_
41#undef FORCE_NOMPI
42}
43namespace rheolef {
44
45// =======================================================================
46// pastix_base_rep
47// =======================================================================
48template<class T, class M>
49class solver_pastix_base_rep : public solver_abstract_rep<T,M> {
50public:
51
52// allocator:
53
54 solver_pastix_base_rep ();
55 explicit solver_pastix_base_rep (const csr<T,M>& a, const solver_option& opt = solver_option());
56 bool initialized() const { return true; }
57 void update_values (const csr<T,M>& a);
58 ~solver_pastix_base_rep ();
59
60// accessors:
61
62 vec<T,M> trans_solve (const vec<T,M>& rhs) const;
63 vec<T,M> solve (const vec<T,M>& rhs) const;
64
65// internal accessors & modifiers:
66protected:
67
68 void load (const csr<T,M>& a, const solver_option& opt = solver_option());
69 bool is_symmetric () const { return _is_sym; }
70 void resize (pastix_int_t n, pastix_int_t nnz);
71 void load_symmetric (const csr<T,M>& a);
72 void load_unsymmetric (const csr<T,M>& a);
73 void load_both_continued (const csr<T,M>& a);
74 void check () const;
75 void symbolic_factorization ();
76 void numeric_factorization ();
77
78// data:
79public: // TODO: protected
80//protected:
81 static const pastix_int_t _base = 1;
82 pastix_int_t _n;
83 pastix_int_t _nnz;
84 mutable std::vector<pastix_int_t> _ptr;
85 mutable std::vector<pastix_int_t> _idx; // row index, in csc format: dis_i = idx[p], p=0..nnz-1
86 mutable std::vector<T> _val;
87 bool _is_sym;
88 size_t _pattern_dimension;
89 mutable pastix_data_t* _pastix_data; // Pointer to a storage structure needed by pastix
90 mutable pastix_int_t _iparm [IPARM_SIZE]; // integer parameters for pastix
91 mutable T _dparm [DPARM_SIZE]; // floating parameters for pastix
92 distributor _csr_row_ownership;
93 distributor _csr_col_ownership;
94 solver_option _opt;
95 mutable std::vector<T> _new_rhs;
96 mutable std::vector<pastix_int_t> _new_i2dis_i_base;
97 mutable std::vector<pastix_int_t> _i2new_dis_i; // permutation
98 bool _have_pastix_bug_small_matrix; // circumvent when np < a.dis_nrow
99 csr<T,M> _a_when_bug; // use it when pastix bug (too small)
100
101// internal:
102 pastix_data_t** pp_data() const { return &_pastix_data; }
103
104};
105// =======================================================================
106// pastix_rep
107// =======================================================================
108template<class T,class M>
109class solver_pastix_rep {};
110
111// ====================================================================
112// sequential pastix interface
113// ====================================================================
114template<class T>
115class solver_pastix_rep<T,sequential> : public solver_pastix_base_rep<T,sequential> {
116public:
117 typedef solver_pastix_base_rep<T,sequential> base;
118
119// allocator:
120
121 solver_pastix_rep () : base() {}
122 explicit solver_pastix_rep (const csr<T,sequential>& a, const solver_option& opt = solver_option())
123 : base(a,opt) {}
124 void update_values (const csr<T,sequential>& a) { base::update_values(a); }
125
126// accessors:
127
128 vec<T,sequential> trans_solve (const vec<T,sequential>& rhs) const { return base::trans_solve(rhs); }
129 vec<T,sequential> solve (const vec<T,sequential>& rhs) const { return base::solve(rhs); }
130
131// internal accessors & modifiers:
132protected:
133
134 void load (const csr<T,sequential>& a, const solver_option& opt = solver_option()) { load(a,opt); }
135 bool is_symmetric () const { return base::is_symmetric(); }
136 void resize (pastix_int_t n, pastix_int_t nnz) { resize (n,nnz); }
137 void load_symmetric (const csr<T,sequential>& a) { load_symmetric(a); }
138 void load_unsymmetric (const csr<T,sequential>& a) { load_unsymmetric(a); }
139 void load_both_continued (const csr<T,sequential>& a) { load_both_continued(a); }
140 void check () const { check(); }
141 void symbolic_factorization () { symbolic_factorization(); }
142 void numeric_factorization () { numeric_factorization(); }
143};
144// ====================================================================
145// distributed pastix interface
146// ====================================================================
147#ifdef _RHEOLEF_HAVE_MPI
148template<class T>
149class solver_pastix_rep<T,distributed> : public solver_pastix_base_rep<T,distributed> {
150 typedef solver_pastix_base_rep<T,distributed> base;
151
152public:
153
154// allocator:
155
156 solver_pastix_rep ();
157 explicit solver_pastix_rep (const csr<T,distributed>& a, const solver_option& opt = solver_option());
158 void update_values (const csr<T,distributed>& a);
159 ~solver_pastix_rep ();
160
161// accessors:
162
163 const communicator& comm () const { return _comm; }
164 vec<T,distributed> trans_solve (const vec<T,distributed>& rhs) const;
165 vec<T,distributed> solve (const vec<T,distributed>& rhs) const;
166
167// internal accessors & modifiers:
168protected:
169
170 void load (const csr<T,distributed>& a, const solver_option& opt = solver_option());
171 bool is_symmetric () const { return base::_is_sym; }
172 void resize (pastix_int_t n, pastix_int_t nnz);
173 void load_symmetric (const csr<T,distributed>& a);
174 void load_unsymmetric (const csr<T,distributed>& a);
175 void load_both_continued (const csr<T,distributed>& a);
176 void check () const;
177 void symbolic_factorization ();
178 void numeric_factorization ();
179
180// data:
181protected:
182 communicator _comm;
183 std::vector<pastix_int_t> _i2dis_i_base; // dis_j = i2dis_i_base[j] - base, j=0..n-1
184 pastix_int_t _new_n; // new re-ordering
185 pastix_int_t* _new_ptr;
186 pastix_int_t* _new_idx;
187 T* _new_val;
188};
189#endif // _RHEOLEF_HAVE_MPI
190
191} // namespace rheolef
192#endif // _RHEOLEF_HAVE_PASTIX
193#endif // _RHEOLEF_SOLVER_PASTIX_H
see the communicator page for the full documentation
distributed
Definition asr.cc:228
Expr1::float_type T
Definition field_expr.h:230
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
void check(Float p, const field &uh)
void load(idiststream &in, Float &p, field &uh)