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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5220 - (show annotations)
Thu Oct 23 03:50:37 2014 UTC (6 years ago) by caltinay
File size: 14259 byte(s)
Fixes to CDS SpMV symmetric and tests for blocksizes 1-4 symm/non-symm.

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

  ViewVC Help
Powered by ViewVC 1.1.26