source: test/unit_test/solver_test.cpp@ 2fad0e0

Last change on this file since 2fad0e0 was 48b662, checked in by Olaf Lenz <olenz@…>, 14 years ago

Moved files in scafacos_fcs one level up.

git-svn-id: ​https://svn.version.fz-juelich.de/scafacos/trunk@847 5161e1c8-67bf-11de-9fd5-51895aff932f
  • Property mode set to 100644
File size: 7.0 KB
Line 
1/*
2 * solver_test.cpp
3 *
4 * Created on: 21.07.2010
5 * Author: Julian Iseringhausen
6 */
7
8#ifdef HAVE_CONFIG_H
9#include <config.h>
10#endif
11
12#include <limits>
13
14#ifdef HAVE_LAPACK
15#include "solver/dgesv.hpp"
16#include "solver/dsysv.hpp"
17#endif
18#include "solver/givens.hpp"
19#include "solver/solver_dirichlet.hpp"
20
21#include "solver_test.hpp"
22
23using namespace VMG;
24using VMGTests::SolverTestSuite;
25
26CPPUNIT_TEST_SUITE_REGISTRATION(SolverTestSuite);
27
28void SolverTestSuite::setUp()
29{
30#ifdef HAVE_LAPACK
31 dgesv = new DGESV<SolverDirichlet>();
32 dsysv = new DSYSV<SolverDirichlet>();
33
34 dgesv->x = new double[5];
35 dsysv->x = new double[10];
36
37 dgesv->vec_size = dgesv->max_size = 5;
38 dsysv->vec_size = dsysv->max_size = 10;
39
40 dgesv->b = new double[5];
41 dsysv->b = new double[10];
42
43 dgesv->A = new double*[5];
44 dsysv->A = new double*[10];
45
46 for (int i=0; i<5; ++i) {
47 dgesv->A[i] = new double[5];
48 }
49
50 for (int i=0; i<10; ++i)
51 dsysv->A[i] = new double[10];
52#endif
53
54 givens = new Givens<SolverDirichlet>();
55 givens->vec_size = givens->max_size = 5;
56 givens->x = new double[5];
57 givens->b = new double[5];
58 givens->A = new double*[5];
59
60 for (int i=0; i<5; ++i) {
61 givens->A[i] = new double[5];
62 }
63}
64
65void SolverTestSuite::tearDown()
66{
67#ifdef HAVE_LAPACK
68 delete dgesv;
69 delete dsysv;
70#endif
71 delete givens;
72}
73
74void SolverTestSuite::DGESVTest()
75{
76#ifdef HAVE_LAPACK
77 dgesv->A[0][0] = 1.0;
78 dgesv->A[0][1] = -4.0;
79 dgesv->A[0][2] = 0.0;
80 dgesv->A[0][3] = -5.0;
81 dgesv->A[0][4] = 0.0;
82
83 dgesv->A[1][0] = 0.0;
84 dgesv->A[1][1] = 2.0;
85 dgesv->A[1][2] = 2.0;
86 dgesv->A[1][3] = 5.0;
87 dgesv->A[1][4] = 4.0;
88
89 dgesv->A[2][0] = 4.0;
90 dgesv->A[2][1] = -16.0;
91 dgesv->A[2][2] = 0.0;
92 dgesv->A[2][3] = -25.0;
93 dgesv->A[2][4] = 0.0;
94
95 dgesv->A[3][0] = -1.0;
96 dgesv->A[3][1] = -2.0;
97 dgesv->A[3][2] = -4.0;
98 dgesv->A[3][3] = -5.0;
99 dgesv->A[3][4] = -8.0;
100
101 dgesv->A[4][0] = 0.0;
102 dgesv->A[4][1] = 0.0;
103 dgesv->A[4][2] = 2.0;
104 dgesv->A[4][3] = 0.0;
105 dgesv->A[4][4] = 5.0;
106
107 dgesv->b[0] = 12.5;
108 dgesv->b[1] = 8.0;
109 dgesv->b[2] = -10.0;
110 dgesv->b[3] = 3.5;
111 dgesv->b[4] = -16.0;
112
113 dgesv->Compute();
114
115 CPPUNIT_ASSERT_DOUBLES_EQUAL( 8.5, dgesv->x[0], 1.0e-12);
116 CPPUNIT_ASSERT_DOUBLES_EQUAL(-16.0, dgesv->x[1], 1.0e-12);
117 CPPUNIT_ASSERT_DOUBLES_EQUAL(-18.0, dgesv->x[2], 1.0e-12);
118 CPPUNIT_ASSERT_DOUBLES_EQUAL( 12.0, dgesv->x[3], 1.0e-12);
119 CPPUNIT_ASSERT_DOUBLES_EQUAL( 4.0, dgesv->x[4], 1.0e-12);
120#endif
121}
122
123
124void SolverTestSuite::DSYSVTest()
125{
126#ifdef HAVE_LAPACK
127 dsysv->A[0][0] = 4.0;
128 dsysv->A[0][1] = -1.0;
129 dsysv->A[0][2] = -1.0;
130 dsysv->A[0][3] = -1.0;
131 dsysv->A[0][4] = 0.0;
132 dsysv->A[0][5] = 0.0;
133 dsysv->A[0][6] = -1.0;
134 dsysv->A[0][7] = 0.0;
135 dsysv->A[0][8] = 0.0;
136 dsysv->A[0][9] = 1.0;
137
138 dsysv->A[1][0] = -1.0;
139 dsysv->A[1][1] = 4.0;
140 dsysv->A[1][2] = -1.0;
141 dsysv->A[1][3] = 0.0;
142 dsysv->A[1][4] = -1.0;
143 dsysv->A[1][5] = 0.0;
144 dsysv->A[1][6] = 0.0;
145 dsysv->A[1][7] = -1.0;
146 dsysv->A[1][8] = 0.0;
147 dsysv->A[1][9] = 1.0;
148
149 dsysv->A[2][0] = -1.0;
150 dsysv->A[2][1] = -1.0;
151 dsysv->A[2][2] = 4.0;
152 dsysv->A[2][3] = 0.0;
153 dsysv->A[2][4] = 0.0;
154 dsysv->A[2][5] = -1.0;
155 dsysv->A[2][6] = 0.0;
156 dsysv->A[2][7] = 0.0;
157 dsysv->A[2][8] = -1.0;
158 dsysv->A[2][9] = 1.0;
159
160 dsysv->A[3][0] = -1.0;
161 dsysv->A[3][1] = 0.0;
162 dsysv->A[3][2] = 0.0;
163 dsysv->A[3][3] = 4.0;
164 dsysv->A[3][4] = -1.0;
165 dsysv->A[3][5] = -1.0;
166 dsysv->A[3][6] = -1.0;
167 dsysv->A[3][7] = 0.0;
168 dsysv->A[3][8] = 0.0;
169 dsysv->A[3][9] = 1.0;
170
171 dsysv->A[4][0] = 0.0;
172 dsysv->A[4][1] = -1.0;
173 dsysv->A[4][2] = 0.0;
174 dsysv->A[4][3] = -1.0;
175 dsysv->A[4][4] = 4.0;
176 dsysv->A[4][5] = -1.0;
177 dsysv->A[4][6] = 0.0;
178 dsysv->A[4][7] = -1.0;
179 dsysv->A[4][8] = 0.0;
180 dsysv->A[4][9] = 1.0;
181
182 dsysv->A[5][0] = 0.0;
183 dsysv->A[5][1] = 0.0;
184 dsysv->A[5][2] = -1.0;
185 dsysv->A[5][3] = -1.0;
186 dsysv->A[5][4] = -1.0;
187 dsysv->A[5][5] = 4.0;
188 dsysv->A[5][6] = 0.0;
189 dsysv->A[5][7] = 0.0;
190 dsysv->A[5][8] = -1.0;
191 dsysv->A[5][9] = 1.0;
192
193 dsysv->A[6][0] = -1.0;
194 dsysv->A[6][1] = 0.0;
195 dsysv->A[6][2] = 0.0;
196 dsysv->A[6][3] = -1.0;
197 dsysv->A[6][4] = 0.0;
198 dsysv->A[6][5] = 0.0;
199 dsysv->A[6][6] = 4.0;
200 dsysv->A[6][7] = -1.0;
201 dsysv->A[6][8] = -1.0;
202 dsysv->A[6][9] = 1.0;
203
204 dsysv->A[7][0] = 0.0;
205 dsysv->A[7][1] = -1.0;
206 dsysv->A[7][2] = 0.0;
207 dsysv->A[7][3] = 0.0;
208 dsysv->A[7][4] = -1.0;
209 dsysv->A[7][5] = 0.0;
210 dsysv->A[7][6] = -1.0;
211 dsysv->A[7][7] = 4.0;
212 dsysv->A[7][8] = -1.0;
213 dsysv->A[7][9] = 1.0;
214
215 dsysv->A[8][0] = 0.0;
216 dsysv->A[8][1] = 0.0;
217 dsysv->A[8][2] = -1.0;
218 dsysv->A[8][3] = 0.0;
219 dsysv->A[8][4] = 0.0;
220 dsysv->A[8][5] = -1.0;
221 dsysv->A[8][6] = -1.0;
222 dsysv->A[8][7] = -1.0;
223 dsysv->A[8][8] = 4.0;
224 dsysv->A[8][9] = 1.0;
225
226 dsysv->A[9][0] = 1.0;
227 dsysv->A[9][1] = 1.0;
228 dsysv->A[9][2] = 1.0;
229 dsysv->A[9][3] = 1.0;
230 dsysv->A[9][4] = 1.0;
231 dsysv->A[9][5] = 1.0;
232 dsysv->A[9][6] = 1.0;
233 dsysv->A[9][7] = 1.0;
234 dsysv->A[9][8] = 1.0;
235 dsysv->A[9][9] = 0.0;
236
237
238 dsysv->b[0] = -0.125;
239 dsysv->b[1] = -0.125;
240 dsysv->b[2] = -0.125;
241 dsysv->b[3] = -0.125;
242 dsysv->b[4] = 1.0;
243 dsysv->b[5] = -0.125;
244 dsysv->b[6] = -0.125;
245 dsysv->b[7] = -0.125;
246 dsysv->b[8] = -0.125;
247 dsysv->b[9] = 0.0;
248
249 dsysv->Compute();
250
251 CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.0625, dsysv->x[0], std::numeric_limits<double>::epsilon());
252 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, dsysv->x[1], std::numeric_limits<double>::epsilon());
253 CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.0625, dsysv->x[2], std::numeric_limits<double>::epsilon());
254 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, dsysv->x[3], std::numeric_limits<double>::epsilon());
255 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.25, dsysv->x[4], std::numeric_limits<double>::epsilon());
256 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, dsysv->x[5], std::numeric_limits<double>::epsilon());
257 CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.0625, dsysv->x[6], std::numeric_limits<double>::epsilon());
258 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, dsysv->x[7], std::numeric_limits<double>::epsilon());
259 CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.0625, dsysv->x[8], std::numeric_limits<double>::epsilon());
260#endif
261}
262
263void SolverTestSuite::GivensTest()
264{
265 givens->A[0][0] = 1.0;
266 givens->A[0][1] = -4.0;
267 givens->A[0][2] = 0.0;
268 givens->A[0][3] = -5.0;
269 givens->A[0][4] = 0.0;
270
271 givens->A[1][0] = 0.0;
272 givens->A[1][1] = 2.0;
273 givens->A[1][2] = 2.0;
274 givens->A[1][3] = 5.0;
275 givens->A[1][4] = 4.0;
276
277 givens->A[2][0] = 4.0;
278 givens->A[2][1] = -16.0;
279 givens->A[2][2] = 0.0;
280 givens->A[2][3] = -25.0;
281 givens->A[2][4] = 0.0;
282
283 givens->A[3][0] = -1.0;
284 givens->A[3][1] = -2.0;
285 givens->A[3][2] = -4.0;
286 givens->A[3][3] = -5.0;
287 givens->A[3][4] = -8.0;
288
289 givens->A[4][0] = 0.0;
290 givens->A[4][1] = 0.0;
291 givens->A[4][2] = 2.0;
292 givens->A[4][3] = 0.0;
293 givens->A[4][4] = 5.0;
294
295 givens->b[0] = 12.5;
296 givens->b[1] = 8.0;
297 givens->b[2] = -10.0;
298 givens->b[3] = 3.5;
299 givens->b[4] = -16.0;
300
301 givens->Compute();
302
303 CPPUNIT_ASSERT_DOUBLES_EQUAL( 8.5, givens->x[0], 1.0e-12);
304 CPPUNIT_ASSERT_DOUBLES_EQUAL(-16.0, givens->x[1], 1.0e-12);
305 CPPUNIT_ASSERT_DOUBLES_EQUAL(-18.0, givens->x[2], 1.0e-12);
306 CPPUNIT_ASSERT_DOUBLES_EQUAL( 12.0, givens->x[3], 1.0e-12);
307 CPPUNIT_ASSERT_DOUBLES_EQUAL( 4.0, givens->x[4], 1.0e-12);
308}
Note: See TracBrowser for help on using the repository browser.