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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (hide 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 caltinay 5212
2     /*****************************************************************************
3     *
4 uqaeller 6939 * Copyright (c) 2003-2020 by The University of Queensland
5 caltinay 5212 * http://www.uq.edu.au
6     *
7     * Primary Business: Queensland, Australia
8 jfenwick 6112 * Licensed under the Apache License, version 2.0
9     * http://www.apache.org/licenses/LICENSE-2.0
10 caltinay 5212 *
11     * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     * Development 2012-2013 by School of Earth Sciences
13 uqaeller 6939 * Development from 2014-2017 by Centre for Geoscience Computing (GeoComp)
14     * Development from 2019 by School of Earth and Environmental Sciences
15     **
16 caltinay 5212 *****************************************************************************/
17    
18     #include "SystemMatrixTestCase.h"
19 caltinay 6001
20 caltinay 5212 #include <ripley/Rectangle.h>
21     #include <ripley/RipleySystemMatrix.h>
22 caltinay 6001
23     #include <escript/FunctionSpaceFactory.h>
24    
25 caltinay 5212 #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 caltinay 5220 /// 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 caltinay 5212 TestSuite* SystemMatrixTestCase::suite()
149     {
150     TestSuite *testSuite = new TestSuite("SystemMatrixTestCase");
151     testSuite->addTest(new TestCaller<SystemMatrixTestCase>(
152 caltinay 5220 "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 caltinay 5212 return testSuite;
176     }
177    
178 caltinay 5220 void SystemMatrixTestCase::setUp()
179 caltinay 5212 {
180 caltinay 6001 mpiInfo = escript::makeInfo(MPI_COMM_WORLD);
181 caltinay 5220 domain.reset(new ripley::Rectangle(4, 3, 0., 0., 1., 1.));
182     }
183 caltinay 5212
184 caltinay 5220 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 caltinay 5212
191 caltinay 5220 // 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 caltinay 5212 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 caltinay 5220 //mat->saveMM("/tmp/test.mtx");
219     return matptr;
220     }
221 caltinay 5212
222 caltinay 5220 escript::Data SystemMatrixTestCase::createInputVector(int blocksize)
223     {
224     escript::FunctionSpace fs(escript::solution(*domain));
225 caltinay 5212 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 caltinay 5220 return x;
235     }
236 caltinay 5212
237 caltinay 5220 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 caltinay 5212
250 caltinay 5220 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 caltinay 5212 }
262    
263 caltinay 5220 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