1 #include "config.h"
2
3 #include "fenv.h"
4 #include "SparseMatrix.h"
5
6 using namespace std;
7 using namespace psurface;
8
9
10 template<typename ctype>
check_result(const Vector<ctype> & result,const Vector<ctype> & trueResult,const ctype & tolerance,char const * const message)11 void check_result(const Vector<ctype>& result,
12 const Vector<ctype>& trueResult,
13 const ctype& tolerance, char const * const message) {
14 const size_t n = trueResult.size();
15 Vector<ctype> res(n);
16
17 if ((result - trueResult).length() > tolerance)
18 throw runtime_error(message);
19 }
20
21 template<typename ctype>
check_result_by_mult(const SparseMatrix<ctype> & matrix,const Vector<ctype> & result,const Vector<ctype> & b,const ctype & tolerance,char const * const message)22 void check_result_by_mult(const SparseMatrix<ctype>& matrix,
23 const Vector<ctype>& result,
24 const Vector<ctype>& b,
25 const ctype& tolerance, char const * const message) {
26 if ((matrix.multVec(result) - b).length() > tolerance)
27 throw runtime_error(message);
28 }
29
30
31 template<typename ctype>
test(const ctype tolerance,const int n,const int maxIter)32 void test(const ctype tolerance, const int n, const int maxIter) {
33 //// setup for n-dimensional tests
34 // id matrix test
35 Vector<ctype> residue(n);
36 Vector<ctype> result(n);
37
38 // setup first rhs
39 Vector<ctype> b(n);
40 for (size_t i = 0; i < n; ++i)
41 b[i] = StaticVector<ctype, 2>(i+1, i+1+n);
42
43 SparseMatrix<ctype> id(n);
44 for (size_t i = 0; i < n; ++i)
45 id.setEntry(i, i, 1);
46
47 id.BiCGSTAB(b, result, residue, maxIter, tolerance);
48 check_result<ctype>(result, b, tolerance, "id matrix check 1 failed");
49
50 // setup another rhs
51 for (size_t i = 0; i < n; ++i)
52 b[i] = StaticVector<ctype, 2>(0, 0);
53
54 id.BiCGSTAB(b, result, residue, maxIter, tolerance);
55 check_result<ctype>(result, b, tolerance, "id matrix check 2 failed");
56
57
58 // ladder matrix test
59 SparseMatrix<ctype> ladder(n);
60 for (size_t i = 0; i < n; ++i)
61 ladder.setEntry(i, i, i+1);
62 for (size_t i = 0; i < n; ++i)
63 b[i] = StaticVector<ctype, 2>(i+1, (i+1)*2);
64
65 // obtain solutions
66 ladder.BiCGSTAB(b, result, residue, maxIter, tolerance);
67 Vector<ctype> trueLadder(n);
68 for (size_t i = 0; i < n; ++i)
69 trueLadder[i] = StaticVector<ctype, 2>(1, 2);
70
71 // check results
72 check_result<ctype>(result, b, tolerance, "ladder matrix check failed");
73
74 ///////////////////////////////////////////////////////////////////////
75 //// setup specific tests
76 ///////////////////////////////////////////////////////////////////////
77
78 Vector<ctype> residue1(3);
79 Vector<ctype> result1(3);
80 Vector<ctype> b1(3);
81 SparseMatrix<ctype> test1(3);
82 /*
83 1 2 3 0 5 -7 -1
84 0 2 0 (x, y) = 4 0 ==> (x, y) = ( 2, 0 )
85 0 0 3 3 6 1 2
86 */
87 test1.setEntry(0, 0, 1);
88 test1.setEntry(0, 1, 2);
89 test1.setEntry(0, 2, 3);
90 test1.setEntry(1, 1, 2);
91 test1.setEntry(2, 2, 3);
92
93 b1[0] = StaticVector<ctype, 2>(0, 5);
94 b1[1] = StaticVector<ctype, 2>(4, 0);
95 b1[2] = StaticVector<ctype, 2>(3, 6);
96
97 // solutions
98 test1.BiCGSTAB(b1, result1, residue1, maxIter, tolerance);
99 Vector<ctype> trueTest1(3);
100 trueTest1[0] = StaticVector<ctype, 2>(-7, -1);
101 trueTest1[1] = StaticVector<ctype, 2>(2, 0);
102 trueTest1[2] = StaticVector<ctype, 2>(1, 2);
103
104 // check results
105 check_result<ctype>(result1, trueTest1, tolerance, "specific matrix check 1 failed");
106
107 // second specific matrix
108 Vector<ctype> residue2(6);
109 Vector<ctype> result2(6);
110 Vector<ctype> b2(6);
111 SparseMatrix<ctype> test2(6);
112
113 test2.setEntry(0, 0, 1);
114 test2.setEntry(1, 1, 1);
115 test2.setEntry(2, 2, 1);
116 test2.setEntry(3, 3, 1);
117 test2.setEntry(4, 4, 1);
118 test2.setEntry(5, 5, 1);
119 test2.setEntry(1, 2, -0.0666667);
120 test2.setEntry(1, 3, -0.312071);
121 test2.setEntry(1, 4, -0.312071);
122 test2.setEntry(1, 0, -0.154595);
123 test2.setEntry(1, 5, -0.154595);
124
125 b2[0] = StaticVector<ctype, 2>(0.5,0.866025);
126 b2[1] = StaticVector<ctype, 2>(0,0);
127 b2[2] = StaticVector<ctype, 2>(1,0);
128 b2[3] = StaticVector<ctype, 2>(-0.5,0.866025);
129 b2[4] = StaticVector<ctype, 2>(-0.5,-0.866025);
130 b2[5] = StaticVector<ctype, 2>(0.5,-0.866025);
131
132 // solutions
133 test2.BiCGSTAB(b2, result2, residue2, maxIter, tolerance);
134
135 // check results
136 check_result_by_mult(test2, result2, b2, tolerance, "specific matrix check 2 failed");
137
138 // third specific matrix -- the identity!
139 Vector<ctype> residual3(2);
140 Vector<ctype> result3(2);
141 Vector<ctype> b3(2);
142 SparseMatrix<ctype> test3(2);
143
144 test3.setEntry(0, 0, 1);
145 test3.setEntry(1, 1, 1);
146
147 b3[0] = StaticVector<ctype, 2>(0,0);
148 b3[1] = StaticVector<ctype, 2>(0,0);
149
150 // solutions
151 test3.BiCGSTAB(b3, result3, residual3, maxIter, tolerance);
152
153 // check results
154 check_result_by_mult(test3, result3, b3, tolerance, "specific matrix check 3 failed");
155 }
156
main(int argc,char * argv[])157 int main (int argc, char* argv[]) {
158
159 feenableexcept(FE_INVALID);
160
161 const int n = 1, maxIter = 3000;
162 const double tolerance = 1e-3;
163
164 try {
165 test<float>(tolerance, n, maxIter);
166 test<double>(tolerance, n, maxIter);
167 test<long double>(tolerance, n, maxIter);
168 } catch (const exception& e) {
169 cout << e.what() << endl;
170
171 return 1;
172 }
173
174 return 0;
175 }
176