/[escript]/release/4.0/ripley/test/SystemMatrixTestCase.cpp
ViewVC logotype

Annotation of /release/4.0/ripley/test/SystemMatrixTestCase.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5220 - (hide annotations)
Thu Oct 23 03:50:37 2014 UTC (6 years, 10 months ago) by caltinay
Original Path: trunk/ripley/test/SystemMatrixTestCase.cpp
File size: 14259 byte(s)
Fixes to CDS SpMV symmetric and tests for blocksizes 1-4 symm/non-symm.

1 caltinay 5212
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 caltinay 5220 /// 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 caltinay 5212 TestSuite* SystemMatrixTestCase::suite()
145     {
146     TestSuite *testSuite = new TestSuite("SystemMatrixTestCase");
147     testSuite->addTest(new TestCaller<SystemMatrixTestCase>(
148 caltinay 5220 "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 caltinay 5212 return testSuite;
172     }
173    
174 caltinay 5220 void SystemMatrixTestCase::setUp()
175 caltinay 5212 {
176 caltinay 5220 mpiInfo = esysUtils::makeInfo(MPI_COMM_WORLD);
177     domain.reset(new ripley::Rectangle(4, 3, 0., 0., 1., 1.));
178     }
179 caltinay 5212
180 caltinay 5220 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 caltinay 5212
187 caltinay 5220 // 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 caltinay 5212 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 caltinay 5220 //mat->saveMM("/tmp/test.mtx");
215     return matptr;
216     }
217 caltinay 5212
218 caltinay 5220 escript::Data SystemMatrixTestCase::createInputVector(int blocksize)
219     {
220     escript::FunctionSpace fs(escript::solution(*domain));
221 caltinay 5212 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 caltinay 5220 return x;
231     }
232 caltinay 5212
233 caltinay 5220 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 caltinay 5212
246 caltinay 5220 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 caltinay 5212 }
258    
259 caltinay 5220 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