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
|
#include <iostream>
#include <vector>
#include <tuple>
#include <amgcl/make_solver.hpp>
#include <amgcl/amg.hpp>
#include <amgcl/backend/mkl.hpp>
#include <amgcl/adapter/crs_builder.hpp>
#include <amgcl/coarsening/plain_aggregates.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/relaxation/spai0.hpp>
#include <amgcl/solver/bicgstab.hpp>
#include <amgcl/profiler.hpp>
#include <tests/sample_problem.hpp>
namespace amgcl {
profiler<> prof;
}
struct poisson_2d {
typedef double val_type;
typedef long col_type;
size_t n;
double h2i;
poisson_2d(size_t n) : n(n), h2i((n - 1) * (n - 1)) {}
size_t rows() const { return n * n; }
size_t nonzeros() const { return 5 * rows(); }
void operator()(size_t row,
std::vector<col_type> &col,
std::vector<val_type> &val
) const
{
size_t i = row % n;
size_t j = row / n;
if (j > 0) {
col.push_back(row - n);
val.push_back(-h2i);
}
if (i > 0) {
col.push_back(row - 1);
val.push_back(-h2i);
}
col.push_back(row);
val.push_back(4 * h2i);
if (i + 1 < n) {
col.push_back(row + 1);
val.push_back(-h2i);
}
if (j + 1 < n) {
col.push_back(row + n);
val.push_back(-h2i);
}
}
};
int main(int argc, char *argv[]) {
using amgcl::prof;
int m = argc > 1 ? atoi(argv[1]) : 1024;
int n = m * m;
#define USE_MKL 1
#if USE_MKL
typedef amgcl::backend::mkl Backend;
amgcl::backend::mkl_vec f(n);
amgcl::backend::mkl_vec x(n);
#else
typedef amgcl::backend::builtin<double> Backend;
std::vector<double> f(n);
std::vector<double> x(n);
#endif
prof.tic("build");
amgcl::make_solver<
amgcl::amg<
Backend,
amgcl::coarsening::smoothed_aggregation,
amgcl::relaxation::spai0
>,
amgcl::solver::bicgstab<Backend>
> solve( amgcl::adapter::make_matrix(poisson_2d(m)) );
prof.toc("build");
//std::cout << solve.amg() << std::endl;
std::fill_n(f.data(), n, 1.0);
std::fill_n(x.data(), n, 0.0);
prof.tic("solve");
size_t iters;
double error;
std::tie(iters, error) = solve(f, x);
prof.toc("solve");
std::cout << "Iterations: " << iters << std::endl
<< "Error: " << error << std::endl
<< std::endl;
std::cout << prof << std::endl;
}
|
