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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
|
#include <iostream>
#include <vector>
#include <map>
#include <amgcl/backend/builtin.hpp>
#include <amgcl/make_solver.hpp>
#include <amgcl/amg.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/coarsening/aggregation.hpp>
#include <amgcl/relaxation/spai0.hpp>
#include <amgcl/solver/cg.hpp>
#include <amgcl/profiler.hpp>
class sparse_matrix {
public:
typedef std::map<int, double> sparse_row;
sparse_matrix(int n, int m) : _n(n), _m(m), _rows(n) { }
int nrows() const { return _n; }
int ncols() const { return _m; }
// Get a value at row i and column j
double operator()(int i, int j) const {
sparse_row::const_iterator elem = _rows[i].find(j);
return elem == _rows[i].end() ? 0.0 : elem->second;
}
// Get reference to a value at row i and column j
double& operator()(int i, int j) { return _rows[i][j]; }
// Access the whole row
const sparse_row& operator[](int i) const { return _rows[i]; }
private:
int _n, _m;
std::vector<sparse_row> _rows;
};
namespace amgcl {
namespace backend {
// Let AMGCL know the value type of our matrix:
template <> struct value_type<sparse_matrix> {
typedef double type;
};
// Let AMGCL know the size of our matrix:
template<> struct rows_impl<sparse_matrix> {
static int get(const sparse_matrix &A) { return A.nrows(); }
};
template<> struct cols_impl<sparse_matrix> {
static int get(const sparse_matrix &A) { return A.ncols(); }
};
template<> struct nonzeros_impl<sparse_matrix> {
static int get(const sparse_matrix &A) {
int n = A.nrows(), nnz = 0;
for(int i = 0; i < n; ++i)
nnz += A[i].size();
return nnz;
}
};
// Allow AMGCL to iterate over the rows of our matrix:
template<> struct row_iterator<sparse_matrix> {
struct iterator {
sparse_matrix::sparse_row::const_iterator _it, _end;
iterator(const sparse_matrix &A, int row)
: _it(A[row].begin()), _end(A[row].end()) { }
// Check if we are at the end of the row.
operator bool() const {
return _it != _end;
}
// Advance to the next nonzero element.
iterator& operator++() {
++_it;
return *this;
}
// Column number of the current nonzero element.
int col() const { return _it->first; }
// Value of the current nonzero element.
double value() const { return _it->second; }
};
typedef iterator type;
};
template<> struct row_begin_impl<sparse_matrix> {
typedef row_iterator<sparse_matrix>::type iterator;
static iterator get(const sparse_matrix &A, int row) {
return iterator(A, row);
}
};
} // namespace backend
profiler<> prof;
} // namespace amgcl
using amgcl::prof;
int main() {
// Discretize a 1D Poisson problem
const int n = 10000;
auto t_total = prof.scoped_tic("total");
sparse_matrix A(n, n);
for(int i = 0; i < n; ++i) {
if (i == 0 || i == n - 1) {
// Dirichlet boundary condition
A(i,i) = 1.0;
} else {
// Internal point.
A(i, i-1) = -1.0;
A(i, i) = 2.0;
A(i, i+1) = -1.0;
}
}
// Create an AMGCL solver for the problem.
typedef amgcl::backend::builtin<double> Backend;
amgcl::make_solver<
amgcl::amg<
Backend,
amgcl::coarsening::aggregation,
amgcl::relaxation::spai0
>,
amgcl::solver::cg<Backend>
> solve( A );
std::cout << solve.precond() << std::endl;
auto t_solve = prof.scoped_tic("solve");
std::vector<double> f(n, 1.0), x(n, 0.0);
solve(f, x);
std::cout << prof << std::endl;
}
|
