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
44{
45namespace ScaledMatrixUnitTests
46{
47typedef NekMatrix<double> InnerMatrix;
48typedef NekMatrix<InnerMatrix, ScaledMatrixTag> SMat;
49
50BOOST_AUTO_TEST_CASE(TestDefaultConstructor)
51{
52 SMat m;
53 BOOST_CHECK_EQUAL(0u, m.GetRows());
54 BOOST_CHECK_EQUAL(0u, m.GetColumns());
55 BOOST_CHECK(std::shared_ptr<InnerMatrix>() != m.GetOwnedMatrix());
56 BOOST_CHECK_EQUAL('N', m.GetTransposeFlag());
57 BOOST_CHECK_EQUAL(0u, m.GetStorageSize());
58 BOOST_CHECK_EQUAL(eFULL, m.GetStorageType());
59 BOOST_CHECK_EQUAL(0.0, m.Scale());
60}
61
62BOOST_AUTO_TEST_CASE(TestConstructorWithInnerMatrix)
63{
64 double buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
65
66 std::shared_ptr<InnerMatrix> in(new InnerMatrix(3, 2, buf));
67 SMat m(2.7, in);
68 BOOST_CHECK_EQUAL(3u, m.GetRows());
69 BOOST_CHECK_EQUAL(2u, m.GetColumns());
70 BOOST_CHECK_EQUAL(in, m.GetOwnedMatrix());
71 BOOST_CHECK_EQUAL('N', m.GetTransposeFlag());
72 BOOST_CHECK_EQUAL(6u, m.GetStorageSize());
73 BOOST_CHECK_EQUAL(eFULL, m.GetStorageType());
74 BOOST_CHECK_EQUAL(2.7, m.Scale());
75}
76
77BOOST_AUTO_TEST_CASE(TestElementAccessLinearAlgerbra)
78{
79 double buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
80
81 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(3, 2, buf));
82 std::shared_ptr<InnerMatrix> in2(new InnerMatrix(3, 2, buf));
83 std::shared_ptr<InnerMatrix> in3(new InnerMatrix(3, 2, buf));
84 std::shared_ptr<InnerMatrix> in4(new InnerMatrix(3, 2, buf));
85 in3->Transpose();
86 in4->Transpose();
87
88 SMat m1(1.0, in1);
89 SMat m2(2.0, in2);
90 SMat m3(3.0, in3);
91 SMat m4(4.0, in4);
92
93 m2.Transpose();
94 m4.Transpose();
95
96 // m1 is N/N
97 // m2 is N/T
98 // m3 is T/N
99 // m4 is T/T
100
101 BOOST_CHECK_EQUAL(m1(0, 0), 1.0);
102 BOOST_CHECK_EQUAL(m1(1, 0), 2.0);
103 BOOST_CHECK_EQUAL(m1(2, 0), 3.0);
104 BOOST_CHECK_EQUAL(m1(0, 1), 4.0);
105 BOOST_CHECK_EQUAL(m1(1, 1), 5.0);
106 BOOST_CHECK_EQUAL(m1(2, 1), 6.0);
107
108 BOOST_CHECK_EQUAL(m4(0, 0), 4.0);
109 BOOST_CHECK_EQUAL(m4(1, 0), 8.0);
110 BOOST_CHECK_EQUAL(m4(2, 0), 12.0);
111 BOOST_CHECK_EQUAL(m4(0, 1), 16.0);
112 BOOST_CHECK_EQUAL(m4(1, 1), 20.0);
113 BOOST_CHECK_EQUAL(m4(2, 1), 24.0);
114
115 BOOST_CHECK_EQUAL(m2(0, 0), 2.0);
116 BOOST_CHECK_EQUAL(m2(0, 1), 4.0);
117 BOOST_CHECK_EQUAL(m2(0, 2), 6.0);
118 BOOST_CHECK_EQUAL(m2(1, 0), 8.0);
119 BOOST_CHECK_EQUAL(m2(1, 1), 10.0);
120 BOOST_CHECK_EQUAL(m2(1, 2), 12.0);
121
122 BOOST_CHECK_EQUAL(m3(0, 0), 3.0);
123 BOOST_CHECK_EQUAL(m3(0, 1), 6.0);
124 BOOST_CHECK_EQUAL(m3(0, 2), 9.0);
125 BOOST_CHECK_EQUAL(m3(1, 0), 12.0);
126 BOOST_CHECK_EQUAL(m3(1, 1), 15.0);
127 BOOST_CHECK_EQUAL(m3(1, 2), 18.0);
128}
129
130BOOST_AUTO_TEST_CASE(TestIteration)
131{
132 double buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
133
134 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(3, 2, buf));
135 std::shared_ptr<InnerMatrix> in2(new InnerMatrix(3, 2, buf));
136 std::shared_ptr<InnerMatrix> in3(new InnerMatrix(3, 2, buf));
137 std::shared_ptr<InnerMatrix> in4(new InnerMatrix(3, 2, buf));
138 in3->Transpose();
139 in4->Transpose();
140
141 SMat m1(1.0, in1);
142 SMat m2(2.0, in2);
143 SMat m3(3.0, in3);
144 SMat m4(4.0, in4);
145
146 m2.Transpose();
147 m4.Transpose();
148
149 // m1 is N/N
150 // m2 is N/T
151 // m3 is T/N
152 // m4 is T/T
153
154 SMat::const_iterator i1 = m1.begin();
155 SMat::const_iterator i2 = m2.begin();
156 SMat::const_iterator i3 = m3.begin();
157 SMat::const_iterator i4 = m4.begin();
158
159 BOOST_CHECK(1.0 == *(i1++));
160 BOOST_CHECK(2.0 == *(i1++));
161 BOOST_CHECK(3.0 == *(i1++));
162 BOOST_CHECK(4.0 == *(i1++));
163 BOOST_CHECK(5.0 == *(i1++));
164 BOOST_CHECK(6.0 == *(i1++));
165 BOOST_CHECK(m1.end() == i1);
166
167 BOOST_CHECK(2.0 == *(i2++));
168 BOOST_CHECK(8.0 == *(i2++));
169 BOOST_CHECK(4.0 == *(i2++));
170 BOOST_CHECK(10.0 == *(i2++));
171 BOOST_CHECK(6.0 == *(i2++));
172 BOOST_CHECK(12.0 == *(i2++));
173 BOOST_CHECK(m2.end() == i2);
174
175 BOOST_CHECK(3.0 == *(i3++));
176 BOOST_CHECK(12.0 == *(i3++));
177 BOOST_CHECK(6.0 == *(i3++));
178 BOOST_CHECK(15.0 == *(i3++));
179 BOOST_CHECK(9.0 == *(i3++));
180 BOOST_CHECK(18.0 == *(i3++));
181 BOOST_CHECK(m3.end() == i3);
182
183 BOOST_CHECK(4.0 == *(i4++));
184 BOOST_CHECK(8.0 == *(i4++));
185 BOOST_CHECK(12.0 == *(i4++));
186 BOOST_CHECK(16.0 == *(i4++));
187 BOOST_CHECK(20.0 == *(i4++));
188 BOOST_CHECK(24.0 == *(i4++));
189 BOOST_CHECK(m4.end() == i4);
190}
191
192BOOST_AUTO_TEST_CASE(TestNMatrixNMatrixMultiplication)
193{
194 {
195 double lhs_buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
196 double rhs_buf[] = {7.0, 8.0, 9.0, 10.0, 11.0, 12.0};
197 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(3, 2, lhs_buf));
198 std::shared_ptr<InnerMatrix> in2(new InnerMatrix(2, 3, rhs_buf));
199
200 SMat m1(1.0, in1);
201 SMat m2(2.0, in2);
202
203 double expected_result_buf[] = {78.0, 108.0, 138.0, 98.0, 136.0,
204 174.0, 118.0, 164.0, 210.0};
205 InnerMatrix expected_result(3, 3, expected_result_buf);
206 InnerMatrix result1 = m1 * m2;
207 BOOST_CHECK_EQUAL(expected_result, result1);
208
209 in1->Transpose();
210 m1.Transpose();
211 InnerMatrix result2 = m1 * m2;
212 BOOST_CHECK_EQUAL(expected_result, result2);
213
214 in2->Transpose();
215 m2.Transpose();
216 InnerMatrix result3 = m1 * m2;
217 BOOST_CHECK_EQUAL(expected_result, result3);
218
219 in1->Transpose();
220 m1.Transpose();
221 InnerMatrix result4 = m1 * m2;
222 BOOST_CHECK_EQUAL(expected_result, result4);
223 }
224}
225
226BOOST_AUTO_TEST_CASE(TestTMatrixNMatrixMultiplication)
227{
228 {
229 double lhs_buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
230 double rhs_buf[] = {7.0, 8.0, 9.0, 10.0, 11.0, 12.0};
231 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(2, 3, lhs_buf));
232 std::shared_ptr<InnerMatrix> in2(new InnerMatrix(2, 3, rhs_buf));
233
234 SMat m1(1.0, in1);
235 SMat m2(2.0, in2);
236
237 double expected_result_buf[] = {46.0, 106.0, 166.0, 58.0, 134.0,
238 210.0, 70.0, 162.0, 254.0};
239 InnerMatrix expected_result(3, 3, expected_result_buf);
240
241 in1->Transpose();
242 InnerMatrix result1 = m1 * m2;
243 BOOST_CHECK_EQUAL(expected_result, result1);
244
245 in1->Transpose();
246 m1.Transpose();
247 InnerMatrix result2 = m1 * m2;
248 BOOST_CHECK_EQUAL(expected_result, result2);
249
250 in1->Transpose();
251 m1.Transpose();
252 in2->Transpose();
253 m2.Transpose();
254 InnerMatrix result3 = m1 * m2;
255 BOOST_CHECK_EQUAL(expected_result, result3);
256
257 in1->Transpose();
258 m1.Transpose();
259 InnerMatrix result4 = m1 * m2;
260 BOOST_CHECK_EQUAL(expected_result, result4);
261 }
262}
263
264BOOST_AUTO_TEST_CASE(TestNMatrixTMatrixMultiplication)
265{
266 {
267 double lhs_buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
268 double rhs_buf[] = {7.0, 8.0, 9.0, 10.0, 11.0, 12.0};
269 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(2, 3, lhs_buf));
270 std::shared_ptr<InnerMatrix> in2(new InnerMatrix(2, 3, rhs_buf));
271
272 SMat m1(1.0, in1);
273 SMat m2(2.0, in2);
274
275 double expected_result_buf[] = {178.0, 232.0, 196.0, 256.0};
276 InnerMatrix expected_result(2, 2, expected_result_buf);
277
278 in2->Transpose();
279 InnerMatrix result1 = m1 * m2;
280 BOOST_CHECK_EQUAL(expected_result, result1);
281
282 in2->Transpose();
283 m2.Transpose();
284 InnerMatrix result2 = m1 * m2;
285 BOOST_CHECK_EQUAL(expected_result, result2);
286
287 in2->Transpose();
288 m2.Transpose();
289 in1->Transpose();
290 m1.Transpose();
291 InnerMatrix result3 = m1 * m2;
292 BOOST_CHECK_EQUAL(expected_result, result3);
293
294 in2->Transpose();
295 m2.Transpose();
296 InnerMatrix result4 = m1 * m2;
297 BOOST_CHECK_EQUAL(expected_result, result4);
298 }
299}
300
301BOOST_AUTO_TEST_CASE(TestTMatrixTMatrixMultiplication)
302{
303 {
304 double lhs_buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
305 double rhs_buf[] = {7.0, 8.0, 9.0, 10.0, 11.0, 12.0};
306 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(3, 2, lhs_buf));
307 std::shared_ptr<InnerMatrix> in2(new InnerMatrix(2, 3, rhs_buf));
308
309 SMat m1(1.0, in1);
310 SMat m2(2.0, in2);
311
312 double expected_result_buf[] = {116.0, 278.0, 128.0, 308.0};
313 InnerMatrix expected_result(2, 2, expected_result_buf);
314
315 in1->Transpose();
316 in2->Transpose();
317 InnerMatrix result1 = m1 * m2;
318 BOOST_CHECK_EQUAL(expected_result, result1);
319
320 in1->Transpose();
321 m1.Transpose();
322 InnerMatrix result2 = m1 * m2;
323 BOOST_CHECK_EQUAL(expected_result, result2);
324
325 in2->Transpose();
326 m2.Transpose();
327 InnerMatrix result3 = m1 * m2;
328 BOOST_CHECK_EQUAL(expected_result, result3);
329
330 in1->Transpose();
331 m1.Transpose();
332 InnerMatrix result4 = m1 * m2;
333 BOOST_CHECK_EQUAL(expected_result, result4);
334 }
335}
336
337BOOST_AUTO_TEST_CASE(TestScaledGlobalTransposeMethod)
338{
339 {
340 double lhs_buf[] = {1, 2, 3, 4, 5, 6};
341
342 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(3, 2, lhs_buf));
343 SMat m1(1, in1);
344 SMat transpose = Transpose(m1);
345
346 BOOST_CHECK_EQUAL(m1(0, 0), transpose(0, 0));
347 BOOST_CHECK_EQUAL(m1(1, 0), transpose(0, 1));
348 BOOST_CHECK_EQUAL(m1(2, 0), transpose(0, 2));
349 BOOST_CHECK_EQUAL(m1(0, 1), transpose(1, 0));
350 BOOST_CHECK_EQUAL(m1(1, 1), transpose(1, 1));
351 BOOST_CHECK_EQUAL(m1(2, 1), transpose(1, 2));
352 }
353}
354
355BOOST_AUTO_TEST_CASE(TestScaledNMatrixVectorMultiply)
356{
357 {
358 double lhs_buf[] = {1.0, 3.0, 5.0, 7.0, 2.0, 4.0, 6.0, 8.0};
359 double rhs_buf[] = {1.0, 2.0};
360
361 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(4, 2, lhs_buf));
362 NekVector<double> rhs(2, rhs_buf);
363 SMat m1(2.0, in1);
364
365 double expected_result_buf[] = {10, 22, 34, 46};
366 NekVector<double> expected_result(4, expected_result_buf);
367
368 NekVector<double> result1 = m1 * rhs;
369 BOOST_CHECK_EQUAL(expected_result, result1);
370
371 in1->Transpose();
372 m1.Transpose();
373 NekVector<double> result2 = m1 * rhs;
374 BOOST_CHECK_EQUAL(expected_result, result2);
375 }
376}
377
378BOOST_AUTO_TEST_CASE(BugReport2)
379{
380 double buf1[] = {1.0, 2.0, 3.0, 4.0};
381 double buf2[] = {5.0, 6.0, 7.0, 8.0};
382 double buf3[] = {-1.0, -2.0, -3.0, -4.0};
383
384 std::shared_ptr<NekMatrix<NekDouble>> DMatInner(
385 new NekMatrix<NekDouble>(2, 2, buf1));
386 std::shared_ptr<NekMatrix<NekDouble>> InvMassInner(
387 new NekMatrix<NekDouble>(2, 2, buf2));
388 NekMatrix<NekDouble> Mat(2, 2, buf3);
389 DNekScalMat DMat(2.0, DMatInner);
390 DNekScalMat InvMass(3.0, InvMassInner);
391
392 Mat = Mat + DMat * InvMass * Transpose(DMat);
393
394 BOOST_CHECK_EQUAL(1391.0, Mat(0, 0));
395 BOOST_CHECK_EQUAL(2037.0, Mat(0, 1));
396 BOOST_CHECK_EQUAL(2062.0, Mat(1, 0));
397 BOOST_CHECK_EQUAL(3020.0, Mat(1, 1));
398}
399
400BOOST_AUTO_TEST_CASE(BugReport1)
401{
402 {
403 double dmat_buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
404 std::shared_ptr<NekMatrix<NekDouble>> inner(
405 new NekMatrix<NekDouble>(2, 3, dmat_buf));
406 DNekScalMat Dmat(2.0, inner);
407
408 std::shared_ptr<NekMatrix<NekDouble>> inner1(
409 new NekMatrix<NekDouble>(2, 3, dmat_buf));
410 DNekScalMat invMass(3.0, inner1);
411
412 NekMatrix<NekDouble> LocMat;
413 LocMat = Transpose(Dmat);
414 LocMat = LocMat * invMass;
415
416 BOOST_CHECK_EQUAL(LocMat(0, 0), 30.0);
417 BOOST_CHECK_EQUAL(LocMat(0, 1), 66.0);
418 BOOST_CHECK_EQUAL(LocMat(0, 2), 102.0);
419 BOOST_CHECK_EQUAL(LocMat(1, 0), 66.0);
420 BOOST_CHECK_EQUAL(LocMat(1, 1), 150.0);
421 BOOST_CHECK_EQUAL(LocMat(1, 2), 234.0);
422 BOOST_CHECK_EQUAL(LocMat(2, 0), 102.0);
423 BOOST_CHECK_EQUAL(LocMat(2, 1), 234.0);
424 BOOST_CHECK_EQUAL(LocMat(2, 2), 366.0);
425 }
426
427 {
428 double dmat_buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
429 std::shared_ptr<NekMatrix<NekDouble>> inner(
430 new NekMatrix<NekDouble>(2, 3, dmat_buf));
431 DNekScalMat Dmat(2.0, inner);
432 std::shared_ptr<NekMatrix<NekDouble>> inner1(
433 new NekMatrix<NekDouble>(2, 3, dmat_buf));
434 DNekScalMat invMass(3.0, inner1);
435
436 NekMatrix<NekDouble> result2 = Transpose(Dmat) * invMass;
437 BOOST_CHECK_EQUAL(result2(0, 0), 30.0);
438 BOOST_CHECK_EQUAL(result2(0, 1), 66.0);
439 BOOST_CHECK_EQUAL(result2(0, 2), 102.0);
440 BOOST_CHECK_EQUAL(result2(1, 0), 66.0);
441 BOOST_CHECK_EQUAL(result2(1, 1), 150.0);
442 BOOST_CHECK_EQUAL(result2(1, 2), 234.0);
443 BOOST_CHECK_EQUAL(result2(2, 0), 102.0);
444 BOOST_CHECK_EQUAL(result2(2, 1), 234.0);
445 BOOST_CHECK_EQUAL(result2(2, 2), 366.0);
446 }
447}
448
449BOOST_AUTO_TEST_CASE(TestScaledTMatrixVectorMultiply)
450{
451 {
452 double lhs_buf[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
453 double rhs_buf[] = {1.0, 2.0};
454
455 std::shared_ptr<InnerMatrix> in1(new InnerMatrix(2, 4, lhs_buf));
456 NekVector<double> rhs(2, rhs_buf);
457 SMat m1(2.0, in1);
458
459 double expected_result_buf[] = {10, 22, 34, 46};
460 NekVector<double> expected_result(4, expected_result_buf);
461
462 SMat m2 = Transpose(m1);
463 NekVector<double> result3 = m2 * rhs;
464 BOOST_CHECK_EQUAL(expected_result, result3);
465
466 in1->Transpose();
467 NekVector<double> result1 = m1 * rhs;
468 BOOST_CHECK_EQUAL(expected_result, result1);
469
470 in1->Transpose();
471 m1.Transpose();
472 NekVector<double> result2 = m1 * rhs;
473 BOOST_CHECK_EQUAL(expected_result, result2);
474 }
475}
476} // namespace ScaledMatrixUnitTests
477} // namespace Nektar
Nektar::ScaledMatrixUnitTests::BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(TestDefaultConstructor)
Definition: TestScaledMatrix.cpp:50
Nektar::ScaledMatrixUnitTests::SMat
NekMatrix< InnerMatrix, ScaledMatrixTag > SMat
Definition: TestScaledMatrix.cpp:48
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:235
Generated on Wed Jan 31 2024 08:26:49 for Nektar++ by doxygen 1.9.4