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
|
#include <cusp/relaxation/jacobi.h>
#include <cusp/relaxation/polynomial.h>
#include <cusp/precond/aggregation/smoothed_aggregation.h>
#include <cusp/krylov/cg.h>
#include <cusp/gallery/poisson.h>
#include <cusp/csr_matrix.h>
#include <iostream>
#include "../timer.h"
template<typename IndexType, typename ValueType, typename MemorySpace>
class unsmoothed_aggregation_options
: public cusp::precond::aggregation::smoothed_aggregation_options<IndexType,ValueType,MemorySpace>
{
typedef cusp::precond::aggregation::smoothed_aggregation_options<IndexType,ValueType,MemorySpace> Parent;
typedef typename Parent::MatrixType MatrixType;
public:
unsmoothed_aggregation_options() : Parent() {}
virtual void smooth_prolongator(const MatrixType& A, const MatrixType& T, MatrixType& P, ValueType& rho_DinvA) const
{
P = T;
}
};
template <typename Monitor>
void report_status(Monitor& monitor)
{
if (monitor.converged())
{
std::cout << "Solver converged to " << monitor.tolerance() << " tolerance";
std::cout << " after " << monitor.iteration_count() << " iterations";
std::cout << " (" << monitor.residual_norm() << " final residual)" << std::endl;
}
else
{
std::cout << "Solver reached iteration limit " << monitor.iteration_limit() << " before converging";
std::cout << " to " << monitor.tolerance() << " tolerance ";
std::cout << " (" << monitor.residual_norm() << " final residual)" << std::endl;
}
}
template<typename MatrixType, typename Prec>
void run_amg(const MatrixType& A, Prec& M)
{
typedef typename MatrixType::index_type IndexType;
typedef typename MatrixType::value_type ValueType;
typedef typename MatrixType::memory_space MemorySpace;
// allocate storage for solution (x) and right hand side (b)
cusp::array1d<ValueType, MemorySpace> x(A.num_rows, 0);
cusp::array1d<ValueType, MemorySpace> b(A.num_rows, 1);
// set stopping criteria (iteration_limit = 1000, relative_tolerance = 1e-10)
cusp::default_monitor<ValueType> monitor(b, 1000, 1e-10);
// solve
timer t1;
cusp::krylov::cg(A, x, b, monitor, M);
std::cout << "solved system in " << t1.milliseconds_elapsed() << " ms " << std::endl;
// report status
report_status(monitor);
}
int main(int argc, char ** argv)
{
typedef int IndexType;
typedef float ValueType;
typedef cusp::device_memory MemorySpace;
// create an empty sparse matrix structure
cusp::hyb_matrix<IndexType, ValueType, MemorySpace> A;
IndexType N = 256;
// create 2D Poisson problem
cusp::gallery::poisson5pt(A, N, N);
// solve without preconditioning
{
std::cout << "\nSolving with no preconditioner" << std::endl;
// allocate storage for solution (x) and right hand side (b)
cusp::array1d<ValueType, MemorySpace> x(A.num_rows, 0);
cusp::array1d<ValueType, MemorySpace> b(A.num_rows, 1);
// set stopping criteria (iteration_limit = 10000, relative_tolerance = 1e-10)
cusp::default_monitor<ValueType> monitor(b, 10000, 1e-10);
// solve
timer t0;
cusp::krylov::cg(A, x, b, monitor);
std::cout << "solved system in " << t0.milliseconds_elapsed() << " ms " << std::endl;
// report status
report_status(monitor);
}
// solve with smoothed aggregation algebraic multigrid preconditioner
{
std::cout << "\nSolving with smoothed aggregation preconditioner and jacobi smoother" << std::endl;
// setup preconditioner
timer t0;
cusp::precond::aggregation::smoothed_aggregation<IndexType, ValueType, MemorySpace> M(A);
std::cout << "constructed hierarchy in " << t0.milliseconds_elapsed() << " ms " << std::endl;
run_amg(A,M);
}
// solve with smoothed aggregation algebraic multigrid preconditioner and polynomial smoother
{
typedef cusp::relaxation::polynomial<ValueType,MemorySpace> Smoother;
std::cout << "\nSolving with smoothed aggregation preconditioner and polynomial smoother" << std::endl;
timer t0;
cusp::precond::aggregation::smoothed_aggregation<IndexType, ValueType, MemorySpace, Smoother> M(A);
std::cout << "constructed hierarchy in " << t0.milliseconds_elapsed() << " ms " << std::endl;
run_amg(A,M);
}
// solve with unsmoothed aggregation algebraic multigrid preconditioner and polynomial smoother
{
std::cout << "\nSolving with unsmoothed aggregation preconditioner" << std::endl;
unsmoothed_aggregation_options<IndexType,ValueType,MemorySpace> opts;
timer t0;
cusp::precond::aggregation::smoothed_aggregation<IndexType, ValueType, MemorySpace> M(A,opts);
std::cout << "constructed hierarchy in " << t0.milliseconds_elapsed() << " ms " << std::endl;
run_amg(A,M);
}
return 0;
}
|
