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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5212 - (show annotations)
Wed Oct 22 05:06:48 2014 UTC (7 years ago) by caltinay
File size: 9050 byte(s)
Preparing ripley SystemMatrix tests.

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 TestSuite* SystemMatrixTestCase::suite()
134 {
135 TestSuite *testSuite = new TestSuite("SystemMatrixTestCase");
136 testSuite->addTest(new TestCaller<SystemMatrixTestCase>(
137 "testSpMV",&SystemMatrixTestCase::testSpMV));
138 return testSuite;
139 }
140
141 void SystemMatrixTestCase::testSpMV()
142 {
143 esysUtils::JMPI info(esysUtils::makeInfo(MPI_COMM_WORLD));
144 escript::Domain_ptr dom(new ripley::Rectangle(4, 3, 0., 0., 1., 1.));
145 escript::FunctionSpace fs(escript::solution(*dom));
146
147 int blocksize = 1;
148 bool symmetric = false;
149 int firstdiag = (symmetric ? 4 : 0);
150 ripley::IndexVector offsets(diag_off+firstdiag, diag_off+9);
151
152 // create a 20x20 matrix with 9 diagonals, blocksize 1, non-symmetric
153 escript::ASM_ptr matptr(new ripley::SystemMatrix(info, blocksize,
154 fs, rows, offsets, symmetric));
155 ripley::SystemMatrix* mat(dynamic_cast<ripley::SystemMatrix*>(matptr.get()));
156
157 ripley::IndexVector rowIdx(4);
158 std::vector<double> array(4*4*blocksize*blocksize);
159
160 for (int i=0; i<8; i++) {
161 rowIdx[0] = 2*i;
162 rowIdx[1] = 2*i+1;
163 rowIdx[2] = 2*i+4;
164 rowIdx[3] = 2*i+5;
165 for (int j=0; j<4*4; j++) {
166 for (int k=0; k<blocksize; k++) {
167 for (int l=0; l<blocksize; l++) {
168 // make main diagonal blocks symmetric since the current
169 // implementation actually reads full main diagonal blocks
170 // so if symmetric flag is set the matrix really has to be
171 // symmetric!
172 array[j*blocksize*blocksize + k*blocksize + l] =
173 (j%5==0 ? 1000*i+50*j+k+l : 1000*i+50*j+blocksize*k+l);
174 }
175 }
176 }
177 mat->add(rowIdx, array);
178 }
179
180 escript::DataTypes::ShapeType shape;
181 if (blocksize > 1)
182 shape.push_back(blocksize);
183 escript::Data x(0., shape, fs, true);
184 for (int i=0; i<rows; i++) {
185 double* xx= x.getSampleDataRW(i);
186 for (int j=0; j<blocksize; j++)
187 xx[j] = (double)i*blocksize + j;
188 }
189
190 //mat->saveMM("/tmp/test.mtx");
191 escript::Data y = mat->vectorMultiply(x);
192 double lsup = 0.;
193 for (int i=0; i<rows; i++) {
194 const double* yy = y.getSampleDataRO(i);
195 const double* yref = (symmetric ? &ref_symm[blocksize-1][i*blocksize] : &ref[blocksize-1][i*blocksize]);
196 for (int j=0; j<blocksize; j++) {
197 lsup = std::max(lsup, std::abs(yref[j] - yy[j]));
198 //std::cerr << yref[j] << " " << yy[j] << std::endl;
199 }
200 }
201
202 CPPUNIT_ASSERT(lsup < 1e-12);
203 }
204

  ViewVC Help
Powered by ViewVC 1.1.26