1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
|
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
///
/// =========================================================================
// preconditioner eigen, seq implementations
//
#include "rheolef/mic.h"
#ifdef _RHEOLEF_HAVE_EIGEN
namespace rheolef {
using namespace std;
// =========================================================================
// the class interface
// =========================================================================
template<class T, class M>
void
solver_mic_rep<T,M>::update_values (const csr<T,M>& a)
{
using namespace Eigen;
Matrix<int,Dynamic,1> nnz_row (a.nrow());
typename csr<T,M>::const_iterator ia = a.begin();
for (size_type i = 0, q = 0, n = a.nrow(); i < n; ++i) {
nnz_row[i] = ia[i+1] - ia[i];
}
SparseMatrix<T> a_tmp (a.nrow(),a.ncol());
if (a.nrow() != 0) {
a_tmp.reserve (nnz_row);
}
for (size_type i = 0, n = a.nrow(); i < n; ++i) {
for (typename csr<T,M>::const_data_iterator p = ia[i]; p < ia[i+1]; ++p) {
a_tmp.insert (i, (*p).first) = (*p).second;
}
}
a_tmp.makeCompressed();
if (a.nrow() != 0) {
_mic_a.setInitialShift (_shift);
_mic_a.compute (a_tmp);
check_macro (_mic_a.info() == Success, "eigen MIC incomplete factorization failed");
}
}
template<class T, class M>
vec<T,M>
solver_mic_rep<T,M>::solve (const vec<T,M>& b) const
{
if (b.dis_size() == 0) return b; // empty matrix
vec<T,M> x(b.ownership());
using namespace Eigen;
Map<Matrix<T,Dynamic,1> > b_map ((T*)(&(*b.begin())), b.size()),
x_map ( &(*x.begin()), x.size());
x_map = _mic_a.solve (b_map);
return x;
}
template<class T, class M>
vec<T,M>
solver_mic_rep<T,M>::trans_solve (const vec<T,M>& b) const
{
return solve(b);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template class solver_mic_rep<Float,sequential>;
#ifdef _RHEOLEF_HAVE_MPI
template class solver_mic_rep<Float,distributed>;
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
#endif // HAVE_EIGEN
|
