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
|
#include <iostream>
#include <vector>
#include <algorithm>
#include <amgcl/amg.hpp>
#include <amgcl/make_solver.hpp>
#include <amgcl/backend/builtin.hpp>
#include <amgcl/adapter/crs_builder.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/relaxation/gauss_seidel.hpp>
#include <amgcl/solver/cg.hpp>
#include <amgcl/profiler.hpp>
#include "sample_problem.hpp"
namespace amgcl {
profiler<> prof;
}
//---------------------------------------------------------------------------
struct poisson_2d {
typedef double val_type;
typedef ptrdiff_t 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);
}
}
};
//---------------------------------------------------------------------------
template <class Vec>
double norm(const Vec &v) {
return sqrt(amgcl::backend::inner_product(v, v));
}
//---------------------------------------------------------------------------
int main(int argc, char *argv[]) {
using amgcl::prof;
int m = argc > 1 ? atoi(argv[1]) : 1024;
int n = m * m;
// Create iterative solver preconditioned by AMG.
// The use of make_matrix() from crs_builder.hpp allows to construct the
// system matrix on demand row by row.
prof.tic("build");
typedef amgcl::make_solver<
amgcl::amg<
amgcl::backend::builtin<double>,
amgcl::coarsening::smoothed_aggregation,
amgcl::relaxation::gauss_seidel
>,
amgcl::solver::cg<
amgcl::backend::builtin<double>
>
> Solver;
Solver solve( amgcl::adapter::make_matrix( poisson_2d(m) ) );
prof.toc("build");
std::cout << solve.precond() << std::endl;
std::vector<double> f(n, 1);
std::vector<double> x(n, 0);
prof.tic("solve");
size_t iters;
double resid;
std::tie(iters, resid) = solve(f, x);
prof.toc("solve");
std::cout << "Solver:" << std::endl
<< " Iterations: " << iters << std::endl
<< " Error: " << resid << std::endl
<< std::endl;
// Use the constructed solver as a preconditioner for another iterative
// solver.
//
// Iterative methods use estimated residual for exit condition. For some
// problems the value of estimated residual can get too far from true
// residual due to round-off errors.
//
// Nesting iterative solvers in this way allows to shave last bits off the
// error.
amgcl::solver::cg< amgcl::backend::builtin<double> > S(n);
std::fill(x.begin(), x.end(), 0);
prof.tic("nested solver");
std::tie(iters, resid) = S(solve.system_matrix(), solve, f, x);
prof.toc("nested solver");
std::cout << "Nested solver:" << std::endl
<< " Iterations: " << iters << std::endl
<< " Error: " << resid << std::endl
<< std::endl;
std::cout << prof << std::endl;
}
|
