/[escript]/trunk/ripley/test/SystemMatrixTestCase.cpp
ViewVC logotype

Contents of /trunk/ripley/test/SystemMatrixTestCase.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (20 months, 3 weeks ago) by uqaeller
File size: 14330 byte(s)
Updated the copyright header.


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2020 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014-2017 by Centre for Geoscience Computing (GeoComp)
14 * Development from 2019 by School of Earth and Environmental Sciences
15 **
16 *****************************************************************************/
17
18 #include "SystemMatrixTestCase.h"
19
20 #include <ripley/Rectangle.h>
21 #include <ripley/RipleySystemMatrix.h>
22
23 #include <escript/FunctionSpaceFactory.h>
24
25 #include <cppunit/TestCaller.h>
26
27 using namespace CppUnit;
28 using namespace std;
29
30 // number of matrix rows (for blocksize 1)
31 const int rows = 20;
32
33 // diagonal offsets for full matrix
34 const int diag_off[] = { -6, -5, -4, -1, 0, 1, 4, 5, 6 };
35
36 // to get these reference values save the matrix via saveMM(), then in python:
37 // print scipy.io.mmread('mat.mtx') * range(20*blocksize)
38
39 // reference results - non-symmetric
40 // block size 1
41 const double ref_bs1[] =
42 { 5800., 6400., 28200., 29200., 72100., 73950.,
43 145600., 148200., 251200., 254600., 388800., 393000.,
44 558400., 563400., 608000., 612800., 360600., 363400.,
45 475800., 479000.};
46 // block size 2
47 const double ref_bs2[] =
48 { 24455., 24507., 27053., 27105., 118083., 118167.,
49 122281., 122365., 299244., 299399., 306990., 307145.,
50 601398., 601614., 612194., 612410., 1031854., 1032134.,
51 1045850., 1046130., 1590310., 1590654., 1607506., 1607850.,
52 2276766., 2277174., 2297162., 2297570., 2475540., 2475931.,
53 2495086., 2495477., 1468199., 1468427., 1479597., 1479825.,
54 1933027., 1933287., 1946025., 1946285.};
55 // block size 3
56 const double ref_bs3[] =
57 { 56174., 56294., 56414., 62156., 62276., 62396.,
58 269954., 270146., 270338., 279536., 279728., 279920.,
59 681940., 682294., 682648., 699604., 699958., 700312.,
60 1368088., 1368580., 1369072., 1392652., 1393144., 1393636.,
61 2342848., 2343484., 2344120., 2374612., 2375248., 2375884.,
62 3605608., 3606388., 3607168., 3644572., 3645352., 3646132.,
63 5156368., 5157292., 5158216., 5202532., 5203456., 5204380.,
64 5603773., 5604658., 5605543., 5647987., 5648872., 5649757.,
65 3323474., 3323990., 3324506., 3349256., 3349772., 3350288.,
66 4372454., 4373042., 4373630., 4401836., 4402424., 4403012.};
67 // block size 4
68 const double ref_bs4[] =
69 { 101334., 101550., 101766., 101982., 112062.,
70 112278., 112494., 112710., 484358., 484702.,
71 485046., 485390., 501486., 501830., 502174.,
72 502518., 1221084., 1221718., 1222352., 1222986.,
73 1252640., 1253274., 1253908., 1254542., 2446892.,
74 2447772., 2448652., 2449532., 2490748., 2491628.,
75 2492508., 2493388., 4185740., 4186876., 4188012.,
76 4189148., 4242396., 4243532., 4244668., 4245804.,
77 6436588., 6437980., 6439372., 6440764., 6506044.,
78 6507436., 6508828., 6510220., 9199436., 9201084.,
79 9202732., 9204380., 9281692., 9283340., 9284988.,
80 9286636., 9994708., 9996286., 9997864., 9999442.,
81 10073464., 10075042., 10076620., 10078198., 5927606.,
82 5928526., 5929446., 5930366., 5973534., 5974454.,
83 5975374., 5976294., 7795430., 7796478., 7797526.,
84 7798574., 7847758., 7848806., 7849854., 7850902.};
85
86 const double* ref[] = {ref_bs1, ref_bs2, ref_bs3, ref_bs4};
87
88 // reference results - symmetric
89 // block size 1
90 const double ref_symm_bs1[] =
91 { 5800., 6400., 28200., 29500., 72100., 75600.,
92 142550., 153300., 243850., 263000., 377150., 404700.,
93 542450., 578400., 587750., 631100., 336050., 382600.,
94 446950., 501200.};
95 // block size 2
96 const double ref_symm_bs2[] =
97 { 24455., 24507., 27202., 27255., 118083., 118167.,
98 123626., 123719., 298942., 299096., 314374., 314574.,
99 588036., 588250., 633214., 633510., 1001284., 1001578.,
100 1080054., 1080446., 1542532., 1542906., 1654894., 1655382.,
101 2211780., 2212234., 2357734., 2358318., 2393346., 2393799.,
102 2568842., 2569441., 1368797., 1369103., 1556820., 1557223.,
103 1816417., 1816771., 2035236., 2035695.};
104 // block size 3
105 const double ref_symm_bs3[] =
106 { 56174., 56294., 56414., 62596., 62722., 62848.,
107 269954., 270146., 270338., 282640., 282874., 283108.,
108 681025., 681376., 681727., 716694., 717246., 717798.,
109 1337096., 1337600., 1338104., 1440210., 1441050., 1441890.,
110 2273084., 2273804., 2274524., 2451726., 2452854., 2453982.,
111 3497072., 3498008., 3498944., 3751242., 3752658., 3754074.,
112 5009060., 5010212., 5011364., 5338758., 5340462., 5342166.,
113 5417693., 5418878., 5420063., 5813769., 5815578., 5817387.,
114 3098622., 3099510., 3100398., 3522842., 3524132., 3525422.,
115 4108830., 4109862., 4110894., 4602314., 4603784., 4605254.};
116 // block size 4
117 const double ref_symm_bs4[] =
118 { 101334., 101550., 101766., 101982., 112920.,
119 113154., 113388., 113622., 484358., 484702.,
120 485046., 485390., 507000., 507458., 507916.,
121 508374., 1219228., 1219856., 1220484., 1221112.,
122 1283176., 1284336., 1285496., 1286656., 2390844.,
123 2391776., 2392708., 2393640., 2575048., 2576848.,
124 2578648., 2580448., 4060604., 4061984., 4063364.,
125 4064744., 4378920., 4381360., 4383800., 4386240.,
126 6242364., 6244192., 6246020., 6247848., 6694792.,
127 6697872., 6700952., 6704032., 8936124., 8938400.,
128 8940676., 8942952., 9522664., 9526384., 9530104.,
129 9533824., 9662308., 9664706., 9667104., 9669502.,
130 10366660., 10370694., 10374728., 10378762., 5526118.,
131 5528050., 5529982., 5531914., 6280848., 6283822.,
132 6286796., 6289770., 7324854., 7327106., 7329358.,
133 7331610., 8202640., 8206030., 8209420., 8212810.};
134
135 const double* ref_symm[] = {ref_symm_bs1, ref_symm_bs2, ref_symm_bs3, ref_symm_bs4};
136
137 /// helper
138 double lsup(const double* d0, const double* d1, int length)
139 {
140 double result = 0.;
141 for (int i=0; i<length; i++) {
142 result = std::max(result, std::abs(d0[i] - d1[i]));
143 //std::cerr << d0[i] << " " << d1[i] << std::endl;
144 }
145 return result;
146 }
147
148 TestSuite* SystemMatrixTestCase::suite()
149 {
150 TestSuite *testSuite = new TestSuite("SystemMatrixTestCase");
151 testSuite->addTest(new TestCaller<SystemMatrixTestCase>(
152 "testSpMV_CPU_blocksize1_nonsymmetric",
153 &SystemMatrixTestCase::testSpMV_CPU_blocksize1_nonsymmetric));
154 testSuite->addTest(new TestCaller<SystemMatrixTestCase>(
155 "testSpMV_CPU_blocksize2_nonsymmetric",
156 &SystemMatrixTestCase::testSpMV_CPU_blocksize2_nonsymmetric));
157 testSuite->addTest(new TestCaller<SystemMatrixTestCase>(
158 "testSpMV_CPU_blocksize3_nonsymmetric",
159 &SystemMatrixTestCase::testSpMV_CPU_blocksize3_nonsymmetric));
160 testSuite->addTest(new TestCaller<SystemMatrixTestCase>(
161 "testSpMV_CPU_blocksize4_nonsymmetric",
162 &SystemMatrixTestCase::testSpMV_CPU_blocksize4_nonsymmetric));
163 testSuite->addTest(new TestCaller<SystemMatrixTestCase>(
164 "testSpMV_CPU_blocksize1_symmetric",
165 &SystemMatrixTestCase::testSpMV_CPU_blocksize1_symmetric));
166 testSuite->addTest(new TestCaller<SystemMatrixTestCase>(
167 "testSpMV_CPU_blocksize2_symmetric",
168 &SystemMatrixTestCase::testSpMV_CPU_blocksize2_symmetric));
169 testSuite->addTest(new TestCaller<SystemMatrixTestCase>(
170 "testSpMV_CPU_blocksize3_symmetric",
171 &SystemMatrixTestCase::testSpMV_CPU_blocksize3_symmetric));
172 testSuite->addTest(new TestCaller<SystemMatrixTestCase>(
173 "testSpMV_CPU_blocksize4_symmetric",
174 &SystemMatrixTestCase::testSpMV_CPU_blocksize4_symmetric));
175 return testSuite;
176 }
177
178 void SystemMatrixTestCase::setUp()
179 {
180 mpiInfo = escript::makeInfo(MPI_COMM_WORLD);
181 domain.reset(new ripley::Rectangle(4, 3, 0., 0., 1., 1.));
182 }
183
184 escript::ASM_ptr SystemMatrixTestCase::createMatrix(int blocksize,
185 bool symmetric)
186 {
187 escript::FunctionSpace fs(escript::solution(*domain));
188 const int firstdiag = (symmetric ? 4 : 0);
189 const ripley::IndexVector offsets(diag_off+firstdiag, diag_off+9);
190
191 // create a matrix with 9 diagonals, given blocksize and symmetric flag
192 escript::ASM_ptr matptr(new ripley::SystemMatrix(mpiInfo, blocksize, fs,
193 rows, offsets, symmetric));
194 ripley::SystemMatrix* mat(dynamic_cast<ripley::SystemMatrix*>(matptr.get()));
195
196 ripley::IndexVector rowIdx(4);
197 std::vector<double> array(4*4*blocksize*blocksize);
198
199 for (int i=0; i<8; i++) {
200 rowIdx[0] = 2*i;
201 rowIdx[1] = 2*i+1;
202 rowIdx[2] = 2*i+4;
203 rowIdx[3] = 2*i+5;
204 for (int j=0; j<4*4; j++) {
205 for (int k=0; k<blocksize; k++) {
206 for (int l=0; l<blocksize; l++) {
207 // make main diagonal blocks symmetric since the current
208 // implementation actually reads full main diagonal blocks
209 // so if symmetric flag is set the matrix really has to be
210 // symmetric!
211 array[j*blocksize*blocksize + k*blocksize + l] =
212 (j%5==0 ? 1000*i+50*j+k+l : 1000*i+50*j+blocksize*k+l);
213 }
214 }
215 }
216 mat->add(rowIdx, array);
217 }
218 //mat->saveMM("/tmp/test.mtx");
219 return matptr;
220 }
221
222 escript::Data SystemMatrixTestCase::createInputVector(int blocksize)
223 {
224 escript::FunctionSpace fs(escript::solution(*domain));
225 escript::DataTypes::ShapeType shape;
226 if (blocksize > 1)
227 shape.push_back(blocksize);
228 escript::Data x(0., shape, fs, true);
229 for (int i=0; i<rows; i++) {
230 double* xx= x.getSampleDataRW(i);
231 for (int j=0; j<blocksize; j++)
232 xx[j] = (double)i*blocksize + j;
233 }
234 return x;
235 }
236
237 void SystemMatrixTestCase::testSpMV_CPU_blocksize1_nonsymmetric()
238 {
239 int blocksize = 1;
240 bool symmetric = false;
241 escript::ASM_ptr mat(createMatrix(blocksize, symmetric));
242 const escript::Data x(createInputVector(blocksize));
243 const escript::Data y(mat->vectorMultiply(x));
244 const double* yref = ref[blocksize-1];
245 const double* yy = y.getSampleDataRO(0);
246 double error = lsup(yref, yy, blocksize*rows);
247 CPPUNIT_ASSERT(error < 1e-12);
248 }
249
250 void SystemMatrixTestCase::testSpMV_CPU_blocksize2_nonsymmetric()
251 {
252 int blocksize = 2;
253 bool symmetric = false;
254 escript::ASM_ptr mat(createMatrix(blocksize, symmetric));
255 const escript::Data x(createInputVector(blocksize));
256 const escript::Data y(mat->vectorMultiply(x));
257 const double* yref = ref[blocksize-1];
258 const double* yy = y.getSampleDataRO(0);
259 double error = lsup(yref, yy, blocksize*rows);
260 CPPUNIT_ASSERT(error < 1e-12);
261 }
262
263 void SystemMatrixTestCase::testSpMV_CPU_blocksize3_nonsymmetric()
264 {
265 int blocksize = 3;
266 bool symmetric = false;
267 escript::ASM_ptr mat(createMatrix(blocksize, symmetric));
268 const escript::Data x(createInputVector(blocksize));
269 const escript::Data y(mat->vectorMultiply(x));
270 const double* yref = ref[blocksize-1];
271 const double* yy = y.getSampleDataRO(0);
272 double error = lsup(yref, yy, blocksize*rows);
273 CPPUNIT_ASSERT(error < 1e-12);
274 }
275
276 void SystemMatrixTestCase::testSpMV_CPU_blocksize4_nonsymmetric()
277 {
278 int blocksize = 4;
279 bool symmetric = false;
280 escript::ASM_ptr mat(createMatrix(blocksize, symmetric));
281 const escript::Data x(createInputVector(blocksize));
282 const escript::Data y(mat->vectorMultiply(x));
283 const double* yref = ref[blocksize-1];
284 const double* yy = y.getSampleDataRO(0);
285 double error = lsup(yref, yy, blocksize*rows);
286 CPPUNIT_ASSERT(error < 1e-12);
287 }
288
289 void SystemMatrixTestCase::testSpMV_CPU_blocksize1_symmetric()
290 {
291 int blocksize = 1;
292 bool symmetric = true;
293 escript::ASM_ptr mat(createMatrix(blocksize, symmetric));
294 const escript::Data x(createInputVector(blocksize));
295 const escript::Data y(mat->vectorMultiply(x));
296 const double* yref = ref_symm[blocksize-1];
297 const double* yy = y.getSampleDataRO(0);
298 double error = lsup(yref, yy, blocksize*rows);
299 CPPUNIT_ASSERT(error < 1e-12);
300 }
301
302 void SystemMatrixTestCase::testSpMV_CPU_blocksize2_symmetric()
303 {
304 int blocksize = 2;
305 bool symmetric = true;
306 escript::ASM_ptr mat(createMatrix(blocksize, symmetric));
307 const escript::Data x(createInputVector(blocksize));
308 const escript::Data y(mat->vectorMultiply(x));
309 const double* yref = ref_symm[blocksize-1];
310 const double* yy = y.getSampleDataRO(0);
311 double error = lsup(yref, yy, blocksize*rows);
312 CPPUNIT_ASSERT(error < 1e-12);
313 }
314
315 void SystemMatrixTestCase::testSpMV_CPU_blocksize3_symmetric()
316 {
317 int blocksize = 3;
318 bool symmetric = true;
319 escript::ASM_ptr mat(createMatrix(blocksize, symmetric));
320 const escript::Data x(createInputVector(blocksize));
321 const escript::Data y(mat->vectorMultiply(x));
322 const double* yref = ref_symm[blocksize-1];
323 const double* yy = y.getSampleDataRO(0);
324 double error = lsup(yref, yy, blocksize*rows);
325 CPPUNIT_ASSERT(error < 1e-12);
326 }
327
328 void SystemMatrixTestCase::testSpMV_CPU_blocksize4_symmetric()
329 {
330 int blocksize = 4;
331 bool symmetric = true;
332 escript::ASM_ptr mat(createMatrix(blocksize, symmetric));
333 const escript::Data x(createInputVector(blocksize));
334 const escript::Data y(mat->vectorMultiply(x));
335 const double* yref = ref_symm[blocksize-1];
336 const double* yy = y.getSampleDataRO(0);
337 double error = lsup(yref, yy, blocksize*rows);
338 CPPUNIT_ASSERT(error < 1e-12);
339 }
340

  ViewVC Help
Powered by ViewVC 1.1.26