Nektar++
TestScaledMatrix.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: TestScaledMatrix.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description:
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <boost/test/unit_test.hpp>
36
37#include <LibUtilities/LinearAlgebra/NekMatrix.hpp>
38#include <boost/test/tools/floating_point_comparison.hpp>
39#include <boost/test/unit_test.hpp>
40
41#include <UnitTests/LibUtilities/LinearAlgebra/TestCombinationRunner.h>
42
43namespace Nektar::ScaledMatrixUnitTests
44{
45typedef NekMatrix<double> InnerMatrix;
46typedef NekMatrix<InnerMatrix, ScaledMatrixTag> SMat;
47
48BOOST_AUTO_TEST_CASE(TestDefaultConstructor)
49{
50 SMat m;
51 BOOST_CHECK_EQUAL(0u, m.GetRows());
52 BOOST_CHECK_EQUAL(0u, m.GetColumns());
53 BOOST_CHECK(std::shared_ptr<InnerMatrix>() != m.GetOwnedMatrix());
54 BOOST_CHECK_EQUAL('N', m.GetTransposeFlag());
55 BOOST_CHECK_EQUAL(0u, m.GetStorageSize());
56 BOOST_CHECK_EQUAL(eFULL, m.GetStorageType());
57 BOOST_CHECK_EQUAL(0.0, m.Scale());
58}
59
60BOOST_AUTO_TEST_CASE(TestConstructorWithInnerMatrix)
61{
62 double buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
63
64 std::shared_ptr<InnerMatrix> in(new InnerMatrix(3, 2, buf));
65 SMat m(2.7, in);
66 BOOST_CHECK_EQUAL(3u, m.GetRows());
67 BOOST_CHECK_EQUAL(2u, m.GetColumns());
68 BOOST_CHECK_EQUAL(in, m.GetOwnedMatrix());
69 BOOST_CHECK_EQUAL('N', m.GetTransposeFlag());
70 BOOST_CHECK_EQUAL(6u, m.GetStorageSize());
71 BOOST_CHECK_EQUAL(eFULL, m.GetStorageType());
72 BOOST_CHECK_EQUAL(2.7, m.Scale());
73}
74
75BOOST_AUTO_TEST_CASE(TestElementAccessLinearAlgerbra)
76{
77 double buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
78
79 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(3, 2, buf));
80 std::shared_ptr<InnerMatrix> in2(new InnerMatrix(3, 2, buf));
81 std::shared_ptr<InnerMatrix> in3(new InnerMatrix(3, 2, buf));
82 std::shared_ptr<InnerMatrix> in4(new InnerMatrix(3, 2, buf));
83 in3->Transpose();
84 in4->Transpose();
85
86 SMat m1(1.0, in1);
87 SMat m2(2.0, in2);
88 SMat m3(3.0, in3);
89 SMat m4(4.0, in4);
90
91 m2.Transpose();
92 m4.Transpose();
93
94 // m1 is N/N
95 // m2 is N/T
96 // m3 is T/N
97 // m4 is T/T
98
99 BOOST_CHECK_EQUAL(m1(0, 0), 1.0);
100 BOOST_CHECK_EQUAL(m1(1, 0), 2.0);
101 BOOST_CHECK_EQUAL(m1(2, 0), 3.0);
102 BOOST_CHECK_EQUAL(m1(0, 1), 4.0);
103 BOOST_CHECK_EQUAL(m1(1, 1), 5.0);
104 BOOST_CHECK_EQUAL(m1(2, 1), 6.0);
105
106 BOOST_CHECK_EQUAL(m4(0, 0), 4.0);
107 BOOST_CHECK_EQUAL(m4(1, 0), 8.0);
108 BOOST_CHECK_EQUAL(m4(2, 0), 12.0);
109 BOOST_CHECK_EQUAL(m4(0, 1), 16.0);
110 BOOST_CHECK_EQUAL(m4(1, 1), 20.0);
111 BOOST_CHECK_EQUAL(m4(2, 1), 24.0);
112
113 BOOST_CHECK_EQUAL(m2(0, 0), 2.0);
114 BOOST_CHECK_EQUAL(m2(0, 1), 4.0);
115 BOOST_CHECK_EQUAL(m2(0, 2), 6.0);
116 BOOST_CHECK_EQUAL(m2(1, 0), 8.0);
117 BOOST_CHECK_EQUAL(m2(1, 1), 10.0);
118 BOOST_CHECK_EQUAL(m2(1, 2), 12.0);
119
120 BOOST_CHECK_EQUAL(m3(0, 0), 3.0);
121 BOOST_CHECK_EQUAL(m3(0, 1), 6.0);
122 BOOST_CHECK_EQUAL(m3(0, 2), 9.0);
123 BOOST_CHECK_EQUAL(m3(1, 0), 12.0);
124 BOOST_CHECK_EQUAL(m3(1, 1), 15.0);
125 BOOST_CHECK_EQUAL(m3(1, 2), 18.0);
126}
127
128BOOST_AUTO_TEST_CASE(TestIteration)
129{
130 double buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
131
132 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(3, 2, buf));
133 std::shared_ptr<InnerMatrix> in2(new InnerMatrix(3, 2, buf));
134 std::shared_ptr<InnerMatrix> in3(new InnerMatrix(3, 2, buf));
135 std::shared_ptr<InnerMatrix> in4(new InnerMatrix(3, 2, buf));
136 in3->Transpose();
137 in4->Transpose();
138
139 SMat m1(1.0, in1);
140 SMat m2(2.0, in2);
141 SMat m3(3.0, in3);
142 SMat m4(4.0, in4);
143
144 m2.Transpose();
145 m4.Transpose();
146
147 // m1 is N/N
148 // m2 is N/T
149 // m3 is T/N
150 // m4 is T/T
151
152 SMat::const_iterator i1 = m1.begin();
153 SMat::const_iterator i2 = m2.begin();
154 SMat::const_iterator i3 = m3.begin();
155 SMat::const_iterator i4 = m4.begin();
156
157 BOOST_CHECK(1.0 == *(i1++));
158 BOOST_CHECK(2.0 == *(i1++));
159 BOOST_CHECK(3.0 == *(i1++));
160 BOOST_CHECK(4.0 == *(i1++));
161 BOOST_CHECK(5.0 == *(i1++));
162 BOOST_CHECK(6.0 == *(i1++));
163 BOOST_CHECK(m1.end() == i1);
164
165 BOOST_CHECK(2.0 == *(i2++));
166 BOOST_CHECK(8.0 == *(i2++));
167 BOOST_CHECK(4.0 == *(i2++));
168 BOOST_CHECK(10.0 == *(i2++));
169 BOOST_CHECK(6.0 == *(i2++));
170 BOOST_CHECK(12.0 == *(i2++));
171 BOOST_CHECK(m2.end() == i2);
172
173 BOOST_CHECK(3.0 == *(i3++));
174 BOOST_CHECK(12.0 == *(i3++));
175 BOOST_CHECK(6.0 == *(i3++));
176 BOOST_CHECK(15.0 == *(i3++));
177 BOOST_CHECK(9.0 == *(i3++));
178 BOOST_CHECK(18.0 == *(i3++));
179 BOOST_CHECK(m3.end() == i3);
180
181 BOOST_CHECK(4.0 == *(i4++));
182 BOOST_CHECK(8.0 == *(i4++));
183 BOOST_CHECK(12.0 == *(i4++));
184 BOOST_CHECK(16.0 == *(i4++));
185 BOOST_CHECK(20.0 == *(i4++));
186 BOOST_CHECK(24.0 == *(i4++));
187 BOOST_CHECK(m4.end() == i4);
188}
189
190BOOST_AUTO_TEST_CASE(TestNMatrixNMatrixMultiplication)
191{
192 {
193 double lhs_buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
194 double rhs_buf[] = {7.0, 8.0, 9.0, 10.0, 11.0, 12.0};
195 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(3, 2, lhs_buf));
196 std::shared_ptr<InnerMatrix> in2(new InnerMatrix(2, 3, rhs_buf));
197
198 SMat m1(1.0, in1);
199 SMat m2(2.0, in2);
200
201 double expected_result_buf[] = {78.0, 108.0, 138.0, 98.0, 136.0,
202 174.0, 118.0, 164.0, 210.0};
203 InnerMatrix expected_result(3, 3, expected_result_buf);
204 InnerMatrix result1 = m1 * m2;
205 BOOST_CHECK_EQUAL(expected_result, result1);
206
207 in1->Transpose();
208 m1.Transpose();
209 InnerMatrix result2 = m1 * m2;
210 BOOST_CHECK_EQUAL(expected_result, result2);
211
212 in2->Transpose();
213 m2.Transpose();
214 InnerMatrix result3 = m1 * m2;
215 BOOST_CHECK_EQUAL(expected_result, result3);
216
217 in1->Transpose();
218 m1.Transpose();
219 InnerMatrix result4 = m1 * m2;
220 BOOST_CHECK_EQUAL(expected_result, result4);
221 }
222}
223
224BOOST_AUTO_TEST_CASE(TestTMatrixNMatrixMultiplication)
225{
226 {
227 double lhs_buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
228 double rhs_buf[] = {7.0, 8.0, 9.0, 10.0, 11.0, 12.0};
229 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(2, 3, lhs_buf));
230 std::shared_ptr<InnerMatrix> in2(new InnerMatrix(2, 3, rhs_buf));
231
232 SMat m1(1.0, in1);
233 SMat m2(2.0, in2);
234
235 double expected_result_buf[] = {46.0, 106.0, 166.0, 58.0, 134.0,
236 210.0, 70.0, 162.0, 254.0};
237 InnerMatrix expected_result(3, 3, expected_result_buf);
238
239 in1->Transpose();
240 InnerMatrix result1 = m1 * m2;
241 BOOST_CHECK_EQUAL(expected_result, result1);
242
243 in1->Transpose();
244 m1.Transpose();
245 InnerMatrix result2 = m1 * m2;
246 BOOST_CHECK_EQUAL(expected_result, result2);
247
248 in1->Transpose();
249 m1.Transpose();
250 in2->Transpose();
251 m2.Transpose();
252 InnerMatrix result3 = m1 * m2;
253 BOOST_CHECK_EQUAL(expected_result, result3);
254
255 in1->Transpose();
256 m1.Transpose();
257 InnerMatrix result4 = m1 * m2;
258 BOOST_CHECK_EQUAL(expected_result, result4);
259 }
260}
261
262BOOST_AUTO_TEST_CASE(TestNMatrixTMatrixMultiplication)
263{
264 {
265 double lhs_buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
266 double rhs_buf[] = {7.0, 8.0, 9.0, 10.0, 11.0, 12.0};
267 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(2, 3, lhs_buf));
268 std::shared_ptr<InnerMatrix> in2(new InnerMatrix(2, 3, rhs_buf));
269
270 SMat m1(1.0, in1);
271 SMat m2(2.0, in2);
272
273 double expected_result_buf[] = {178.0, 232.0, 196.0, 256.0};
274 InnerMatrix expected_result(2, 2, expected_result_buf);
275
276 in2->Transpose();
277 InnerMatrix result1 = m1 * m2;
278 BOOST_CHECK_EQUAL(expected_result, result1);
279
280 in2->Transpose();
281 m2.Transpose();
282 InnerMatrix result2 = m1 * m2;
283 BOOST_CHECK_EQUAL(expected_result, result2);
284
285 in2->Transpose();
286 m2.Transpose();
287 in1->Transpose();
288 m1.Transpose();
289 InnerMatrix result3 = m1 * m2;
290 BOOST_CHECK_EQUAL(expected_result, result3);
291
292 in2->Transpose();
293 m2.Transpose();
294 InnerMatrix result4 = m1 * m2;
295 BOOST_CHECK_EQUAL(expected_result, result4);
296 }
297}
298
299BOOST_AUTO_TEST_CASE(TestTMatrixTMatrixMultiplication)
300{
301 {
302 double lhs_buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
303 double rhs_buf[] = {7.0, 8.0, 9.0, 10.0, 11.0, 12.0};
304 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(3, 2, lhs_buf));
305 std::shared_ptr<InnerMatrix> in2(new InnerMatrix(2, 3, rhs_buf));
306
307 SMat m1(1.0, in1);
308 SMat m2(2.0, in2);
309
310 double expected_result_buf[] = {116.0, 278.0, 128.0, 308.0};
311 InnerMatrix expected_result(2, 2, expected_result_buf);
312
313 in1->Transpose();
314 in2->Transpose();
315 InnerMatrix result1 = m1 * m2;
316 BOOST_CHECK_EQUAL(expected_result, result1);
317
318 in1->Transpose();
319 m1.Transpose();
320 InnerMatrix result2 = m1 * m2;
321 BOOST_CHECK_EQUAL(expected_result, result2);
322
323 in2->Transpose();
324 m2.Transpose();
325 InnerMatrix result3 = m1 * m2;
326 BOOST_CHECK_EQUAL(expected_result, result3);
327
328 in1->Transpose();
329 m1.Transpose();
330 InnerMatrix result4 = m1 * m2;
331 BOOST_CHECK_EQUAL(expected_result, result4);
332 }
333}
334
335BOOST_AUTO_TEST_CASE(TestScaledGlobalTransposeMethod)
336{
337 {
338 double lhs_buf[] = {1, 2, 3, 4, 5, 6};
339
340 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(3, 2, lhs_buf));
341 SMat m1(1, in1);
342 SMat transpose = Transpose(m1);
343
344 BOOST_CHECK_EQUAL(m1(0, 0), transpose(0, 0));
345 BOOST_CHECK_EQUAL(m1(1, 0), transpose(0, 1));
346 BOOST_CHECK_EQUAL(m1(2, 0), transpose(0, 2));
347 BOOST_CHECK_EQUAL(m1(0, 1), transpose(1, 0));
348 BOOST_CHECK_EQUAL(m1(1, 1), transpose(1, 1));
349 BOOST_CHECK_EQUAL(m1(2, 1), transpose(1, 2));
350 }
351}
352
353BOOST_AUTO_TEST_CASE(TestScaledNMatrixVectorMultiply)
354{
355 {
356 double lhs_buf[] = {1.0, 3.0, 5.0, 7.0, 2.0, 4.0, 6.0, 8.0};
357 double rhs_buf[] = {1.0, 2.0};
358
359 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(4, 2, lhs_buf));
360 NekVector<double> rhs(2, rhs_buf);
361 SMat m1(2.0, in1);
362
363 double expected_result_buf[] = {10, 22, 34, 46};
364 NekVector<double> expected_result(4, expected_result_buf);
365
366 NekVector<double> result1 = m1 * rhs;
367 BOOST_CHECK_EQUAL(expected_result, result1);
368
369 in1->Transpose();
370 m1.Transpose();
371 NekVector<double> result2 = m1 * rhs;
372 BOOST_CHECK_EQUAL(expected_result, result2);
373 }
374}
375
376BOOST_AUTO_TEST_CASE(BugReport2)
377{
378 double buf1[] = {1.0, 2.0, 3.0, 4.0};
379 double buf2[] = {5.0, 6.0, 7.0, 8.0};
380 double buf3[] = {-1.0, -2.0, -3.0, -4.0};
381
382 std::shared_ptr<NekMatrix<NekDouble>> DMatInner(
383 new NekMatrix<NekDouble>(2, 2, buf1));
384 std::shared_ptr<NekMatrix<NekDouble>> InvMassInner(
385 new NekMatrix<NekDouble>(2, 2, buf2));
386 NekMatrix<NekDouble> Mat(2, 2, buf3);
387 DNekScalMat DMat(2.0, DMatInner);
388 DNekScalMat InvMass(3.0, InvMassInner);
389
390 Mat = Mat + DMat * InvMass * Transpose(DMat);
391
392 BOOST_CHECK_EQUAL(1391.0, Mat(0, 0));
393 BOOST_CHECK_EQUAL(2037.0, Mat(0, 1));
394 BOOST_CHECK_EQUAL(2062.0, Mat(1, 0));
395 BOOST_CHECK_EQUAL(3020.0, Mat(1, 1));
396}
397
398BOOST_AUTO_TEST_CASE(BugReport1)
399{
400 {
401 double dmat_buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
402 std::shared_ptr<NekMatrix<NekDouble>> inner(
403 new NekMatrix<NekDouble>(2, 3, dmat_buf));
404 DNekScalMat Dmat(2.0, inner);
405
406 std::shared_ptr<NekMatrix<NekDouble>> inner1(
407 new NekMatrix<NekDouble>(2, 3, dmat_buf));
408 DNekScalMat invMass(3.0, inner1);
409
410 NekMatrix<NekDouble> LocMat;
411 LocMat = Transpose(Dmat);
412 LocMat = LocMat * invMass;
413
414 BOOST_CHECK_EQUAL(LocMat(0, 0), 30.0);
415 BOOST_CHECK_EQUAL(LocMat(0, 1), 66.0);
416 BOOST_CHECK_EQUAL(LocMat(0, 2), 102.0);
417 BOOST_CHECK_EQUAL(LocMat(1, 0), 66.0);
418 BOOST_CHECK_EQUAL(LocMat(1, 1), 150.0);
419 BOOST_CHECK_EQUAL(LocMat(1, 2), 234.0);
420 BOOST_CHECK_EQUAL(LocMat(2, 0), 102.0);
421 BOOST_CHECK_EQUAL(LocMat(2, 1), 234.0);
422 BOOST_CHECK_EQUAL(LocMat(2, 2), 366.0);
423 }
424
425 {
426 double dmat_buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
427 std::shared_ptr<NekMatrix<NekDouble>> inner(
428 new NekMatrix<NekDouble>(2, 3, dmat_buf));
429 DNekScalMat Dmat(2.0, inner);
430 std::shared_ptr<NekMatrix<NekDouble>> inner1(
431 new NekMatrix<NekDouble>(2, 3, dmat_buf));
432 DNekScalMat invMass(3.0, inner1);
433
434 NekMatrix<NekDouble> result2 = Transpose(Dmat) * invMass;
435 BOOST_CHECK_EQUAL(result2(0, 0), 30.0);
436 BOOST_CHECK_EQUAL(result2(0, 1), 66.0);
437 BOOST_CHECK_EQUAL(result2(0, 2), 102.0);
438 BOOST_CHECK_EQUAL(result2(1, 0), 66.0);
439 BOOST_CHECK_EQUAL(result2(1, 1), 150.0);
440 BOOST_CHECK_EQUAL(result2(1, 2), 234.0);
441 BOOST_CHECK_EQUAL(result2(2, 0), 102.0);
442 BOOST_CHECK_EQUAL(result2(2, 1), 234.0);
443 BOOST_CHECK_EQUAL(result2(2, 2), 366.0);
444 }
445}
446
447BOOST_AUTO_TEST_CASE(TestScaledTMatrixVectorMultiply)
448{
449 {
450 double lhs_buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
451 double rhs_buf[] = {1.0, 2.0};
452
453 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(2, 4, lhs_buf));
454 NekVector<double> rhs(2, rhs_buf);
455 SMat m1(2.0, in1);
456
457 double expected_result_buf[] = {10, 22, 34, 46};
458 NekVector<double> expected_result(4, expected_result_buf);
459
460 SMat m2 = Transpose(m1);
461 NekVector<double> result3 = m2 * rhs;
462 BOOST_CHECK_EQUAL(expected_result, result3);
463
464 in1->Transpose();
465 NekVector<double> result1 = m1 * rhs;
466 BOOST_CHECK_EQUAL(expected_result, result1);
467
468 in1->Transpose();
469 m1.Transpose();
470 NekVector<double> result2 = m1 * rhs;
471 BOOST_CHECK_EQUAL(expected_result, result2);
472 }
473}
474} // namespace Nektar::ScaledMatrixUnitTests
Nektar::ScaledMatrixUnitTests::BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(TestDefaultConstructor)
Definition: TestScaledMatrix.cpp:48
Nektar::ScaledMatrixUnitTests::SMat
NekMatrix< InnerMatrix, ScaledMatrixTag > SMat
Definition: TestScaledMatrix.cpp:46
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:235
Generated on Wed Oct 30 2024 16:07:32 for Nektar++ by doxygen 1.9.4