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
|
program solve;
{$APPTYPE CONSOLE}
uses
SysUtils,
amgcl in 'amgcl.pas';
const
m = 256;
n = m * m;
nnz = m * m + 4 * (m - 2) * (m - 2);
var
i, j, head, idx: Integer;
ptr, col: Array of Integer;
val, rhs, x: Array of Double;
prm: amgcl.TParams;
solver: amgcl.TSolver;
conv: amgcl.TConvInfo;
begin
try
amgcl.load;
except
on e: Exception do begin
writeln(e.Message);
exit;
end;
end;
SetLength(ptr, n + 1);
SetLength(col, nnz);
SetLength(val, nnz);
SetLength(rhs, n);
SetLength(x, n);
ptr[0] := 0;
idx := 0;
head := 0;
for j := 0 to m - 1 do begin
for i := 0 to m - 1 do begin
if (i = 0) or (i = m - 1) or (j = 0) or (j = m - 1) then
begin
col[head] := idx;
val[head] := 1;
rhs[idx] := 0;
head := head + 1;
end else begin
col[head+0] := idx - m;
col[head+1] := idx - 1;
col[head+2] := idx;
col[head+3] := idx + 1;
col[head+4] := idx + m;
val[head+0] := -1;
val[head+1] := -1;
val[head+2] := 4;
val[head+3] := -1;
val[head+4] := -1;
rhs[idx] := 1;
head := head + 5;
end;
x[idx] := 0;
idx := idx + 1;
ptr[idx] := head;
end
end;
prm := amgcl.TParams.Create;
solver := amgcl.TSolver.Create(
amgcl.coarseningSmoothedAggregation,
amgcl.relaxationSPAI0,
amgcl.solverBiCGStabL,
prm, n, ptr, col, val
);
conv := solver.solve(rhs, x);
writeln('iters: ', conv.iterations);
writeln('resid: ', conv.residual);
solver.free;
prm.free;
amgcl.unload;
end.
|
