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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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     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