Nektar++
Loading...
Searching...
No Matches
testNekMatrix.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: testNekMatrix.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: Tests NekMatrix functionality.
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <LibUtilities/LinearAlgebra/NekMatrix.hpp>
36#include <UnitTests/util.h>
37#include <boost/test/tools/floating_point_comparison.hpp>
38#include <boost/test/unit_test.hpp>
39#include <iostream>
40
41namespace Nektar
42{
43namespace MatrixUnitTests
44{
45BOOST_AUTO_TEST_CASE(TestFullMatrixInversion)
46{
47 {
48 double buf[] = {1.0, 3.0, 2.0, 4.0};
49 NekMatrix<double> m(2, 2, buf);
50 m.Invert();
51
52 BOOST_CHECK_CLOSE(*m(0, 0), -2.0, .00001);
53 BOOST_CHECK_CLOSE(*m(0, 1), 1.0, .00001);
54 BOOST_CHECK_CLOSE(*m(1, 0), 3.0 / 2.0, .00001);
55 BOOST_CHECK_CLOSE(*m(1, 1), -1.0 / 2.0, .00001);
56 }
57
58 {
59 double buf[] = {1.7, 5.6, 9.1, -3.4, -2.0, 4.5, 7.8, 8.2, -2.5,
60 5.0, 9.0, 3.6, 7.3, .98, 1.0, -12.6, 2.9, 6.4,
61 .02, -4.0, -1, 1.7, 5.6, 2.98, 5.0};
62 NekMatrix<double> m(5, 5, buf);
63 m.Invert();
64
65 BOOST_CHECK_CLOSE(*m(0, 0), 0.0005010995978, .000001);
66 BOOST_CHECK_CLOSE(*m(0, 1), -0.4704403712, .000001);
67 BOOST_CHECK_CLOSE(*m(0, 2), 0.2719063614, .000001);
68 BOOST_CHECK_CLOSE(*m(0, 3), -0.3941557805, .000001);
69 BOOST_CHECK_CLOSE(*m(0, 4), 0.09043166650, .000001);
70
71 BOOST_CHECK_CLOSE(*m(1, 0), -0.01087166322, .000001);
72 BOOST_CHECK_CLOSE(*m(1, 1), 0.3242048735, .000001);
73 BOOST_CHECK_CLOSE(*m(1, 2), -0.1605116333, .000001);
74 BOOST_CHECK_CLOSE(*m(1, 3), 0.09133673974, .000001);
75 BOOST_CHECK_CLOSE(*m(1, 4), 0.01293234281, .000001);
76
77 BOOST_CHECK_CLOSE(*m(2, 0), 0.06465598241, .000001);
78 BOOST_CHECK_CLOSE(*m(2, 1), 0.2593895742, .000001);
79 BOOST_CHECK_CLOSE(*m(2, 2), -0.09057233795, .000001);
80 BOOST_CHECK_CLOSE(*m(2, 3), 0.3106562459, .000001);
81 BOOST_CHECK_CLOSE(*m(2, 4), -0.1589713628, .000001);
82
83 BOOST_CHECK_CLOSE(*m(3, 0), -0.03464982687, .000001);
84 BOOST_CHECK_CLOSE(*m(3, 1), 0.2655177914, .000001);
85 BOOST_CHECK_CLOSE(*m(3, 2), -0.1016866584, .000001);
86 BOOST_CHECK_CLOSE(*m(3, 3), 0.2125363445, .000001);
87 BOOST_CHECK_CLOSE(*m(3, 4), -0.1099886183, .000001);
88
89 BOOST_CHECK_CLOSE(*m(4, 0), -0.02957895492, .000001);
90 BOOST_CHECK_CLOSE(*m(4, 1), -0.3518447037, .000001);
91 BOOST_CHECK_CLOSE(*m(4, 2), 0.2060393188, .000001);
92 BOOST_CHECK_CLOSE(*m(4, 3), -0.1411012255, .000001);
93 BOOST_CHECK_CLOSE(*m(4, 4), 0.1670437017, .000001);
94 }
95
96 {
97 }
98}
99
100BOOST_AUTO_TEST_CASE(TestFullMatrixPseudoInversion)
101{
102 double buf[] = {1.0, 3.0, 5.0, 2.0, 4.0, 6.0};
103 NekMatrix<double> m(3, 2, buf);
104
105 m.Invert();
106
107 BOOST_CHECK_EQUAL(m.GetRows(), 2u);
108 BOOST_CHECK_EQUAL(m.GetColumns(), 3u);
109
110 double epsilon = 1.0e-10;
111 BOOST_CHECK_CLOSE(*m(0, 0), -4.0 / 3.0, epsilon);
112 BOOST_CHECK_CLOSE(*m(0, 1), -1.0 / 3.0, epsilon);
113 BOOST_CHECK_CLOSE(*m(0, 2), 2.0 / 3.0, epsilon);
114
115 BOOST_CHECK_CLOSE(*m(1, 0), 13.0 / 12.0, epsilon);
116 BOOST_CHECK_CLOSE(*m(1, 1), 1.0 / 3.0, epsilon);
117 BOOST_CHECK_CLOSE(*m(1, 2), -5.0 / 12.0, epsilon);
118}
119
120BOOST_AUTO_TEST_CASE(TestDiagonalMatrixInversion)
121{
122 double buf[] = {1.0, 2.0, 3.0, 4.0};
123 NekMatrix<double> m(4, 4, buf, eDIAGONAL);
124 m.Invert();
125
126 BOOST_CHECK_EQUAL(m(0, 0), 1.0 / 1.0);
127 BOOST_CHECK_EQUAL(m(1, 1), 1.0 / 2.0);
128 BOOST_CHECK_EQUAL(m(2, 2), 1.0 / 3.0);
129 BOOST_CHECK_EQUAL(m(3, 3), 1.0 / 4.0);
130}
131
132BOOST_AUTO_TEST_CASE(TestNekMatrixConstruction)
133{
134 {
135 std::shared_ptr<Nektar::Matrix<double>> a(
136 new Nektar::NekMatrix<double>(3, 4));
137 std::shared_ptr<Nektar::NekMatrix<double>> b(
138 new Nektar::NekMatrix<double>(5, 6));
139
140 BOOST_CHECK_EQUAL(a->GetRows(), 3u);
141 BOOST_CHECK_EQUAL(a->GetColumns(), 4u);
142 BOOST_CHECK_EQUAL(a->GetSize()[0], 3u);
143 BOOST_CHECK_EQUAL(a->GetSize()[1], 4u);
144
145 BOOST_CHECK_EQUAL(b->GetRows(), 5u);
146 BOOST_CHECK_EQUAL(b->GetColumns(), 6u);
147 BOOST_CHECK_EQUAL(b->GetSize()[0], 5u);
148 BOOST_CHECK_EQUAL(b->GetSize()[1], 6u);
149
150 BOOST_CHECK_EQUAL(a->GetStorageType(), eFULL);
151 BOOST_CHECK_EQUAL(b->GetStorageType(), eFULL);
152
153 BOOST_CHECK_EQUAL(a->GetStorageSize(), 3 * 4);
154 BOOST_CHECK_EQUAL(b->GetStorageSize(), 5 * 6);
155 }
156
157 {
158 std::shared_ptr<Nektar::Matrix<double>> a(
159 new Nektar::NekMatrix<double>(3, 3, eDIAGONAL));
160 std::shared_ptr<Nektar::NekMatrix<double>> b(
161 new Nektar::NekMatrix<double>(5, 5, eDIAGONAL));
162
163 BOOST_CHECK_EQUAL(a->GetRows(), 3);
164 BOOST_CHECK_EQUAL(a->GetColumns(), 3);
165 BOOST_CHECK_EQUAL(a->GetSize()[0], 3);
166 BOOST_CHECK_EQUAL(a->GetSize()[1], 3);
167
168 BOOST_CHECK_EQUAL(b->GetRows(), 5);
169 BOOST_CHECK_EQUAL(b->GetColumns(), 5);
170 BOOST_CHECK_EQUAL(b->GetSize()[0], 5);
171 BOOST_CHECK_EQUAL(b->GetSize()[1], 5);
172
173 BOOST_CHECK_EQUAL(a->GetStorageType(), eDIAGONAL);
174 BOOST_CHECK_EQUAL(b->GetStorageType(), eDIAGONAL);
175
176 BOOST_CHECK_EQUAL(a->GetStorageSize(), 3);
177 BOOST_CHECK_EQUAL(b->GetStorageSize(), 5);
178 }
179}
180
181BOOST_AUTO_TEST_CASE(TestFullNekMatrixGetValue)
182{
183 UnitTests::RedirectCerrIfNeeded();
184 double data[] = {1.0, 89.0, 0.0, 45.12, 2.0, -12.3,
185 892.2532, 76.12, 3.0, -56.7, 211.22, 45.23};
186
187 Nektar::NekMatrix<double> m(4, 3, data);
188
189 BOOST_CHECK_EQUAL(m(0, 0), 1.0);
190 BOOST_CHECK_EQUAL(m(0, 1), 2.0);
191 BOOST_CHECK_EQUAL(m(0, 2), 3.0);
192
193 BOOST_CHECK_EQUAL(m(1, 0), 89.0);
194 BOOST_CHECK_EQUAL(m(1, 1), -12.3);
195 BOOST_CHECK_EQUAL(m(1, 2), -56.7);
196
197 BOOST_CHECK_EQUAL(m(2, 0), 0.0);
198 BOOST_CHECK_EQUAL(m(2, 1), 892.2532);
199 BOOST_CHECK_EQUAL(m(2, 2), 211.22);
200
201 BOOST_CHECK_EQUAL(m(3, 0), 45.12);
202 BOOST_CHECK_EQUAL(m(3, 1), 76.12);
203 BOOST_CHECK_EQUAL(m(3, 2), 45.23);
204
205#ifdef NEKTAR_FULLDEBUG
206 BOOST_CHECK_THROW(m(0, 3), ErrorUtil::NekError);
207 BOOST_CHECK_THROW(m(4, 1), ErrorUtil::NekError);
208 BOOST_CHECK_THROW(m(4, 3), ErrorUtil::NekError);
209#endif // NEKTAR_FULLDEBUG
210
211 std::shared_ptr<Nektar::Matrix<double>> m1(
212 new Nektar::NekMatrix<double>(4, 3, data));
213 std::shared_ptr<Nektar::NekMatrix<double>> m2(
214 new Nektar::NekMatrix<double>(4, 3, data));
215
216 BOOST_CHECK_EQUAL((*m1)(0, 0), 1.0);
217 BOOST_CHECK_EQUAL((*m1)(0, 1), 2.0);
218 BOOST_CHECK_EQUAL((*m1)(0, 2), 3.0);
219
220 BOOST_CHECK_EQUAL((*m1)(1, 0), 89.0);
221 BOOST_CHECK_EQUAL((*m1)(1, 1), -12.3);
222 BOOST_CHECK_EQUAL((*m1)(1, 2), -56.7);
223
224 BOOST_CHECK_EQUAL((*m1)(2, 0), 0.0);
225 BOOST_CHECK_EQUAL((*m1)(2, 1), 892.2532);
226 BOOST_CHECK_EQUAL((*m1)(2, 2), 211.22);
227
228 BOOST_CHECK_EQUAL((*m1)(3, 0), 45.12);
229 BOOST_CHECK_EQUAL((*m1)(3, 1), 76.12);
230 BOOST_CHECK_EQUAL((*m1)(3, 2), 45.23);
231
232#ifdef NEKTAR_FULLDEBUG
233 BOOST_CHECK_THROW((*m1)(0, 3), ErrorUtil::NekError);
234 BOOST_CHECK_THROW((*m1)(4, 1), ErrorUtil::NekError);
235 BOOST_CHECK_THROW((*m1)(4, 3), ErrorUtil::NekError);
236#endif // NEKTAR_FULLDEBUG
237
238 BOOST_CHECK_EQUAL((*m2)(0, 0), 1.0);
239 BOOST_CHECK_EQUAL((*m2)(0, 1), 2.0);
240 BOOST_CHECK_EQUAL((*m2)(0, 2), 3.0);
241
242 BOOST_CHECK_EQUAL((*m2)(1, 0), 89.0);
243 BOOST_CHECK_EQUAL((*m2)(1, 1), -12.3);
244 BOOST_CHECK_EQUAL((*m2)(1, 2), -56.7);
245
246 BOOST_CHECK_EQUAL((*m2)(2, 0), 0.0);
247 BOOST_CHECK_EQUAL((*m2)(2, 1), 892.2532);
248 BOOST_CHECK_EQUAL((*m2)(2, 2), 211.22);
249
250 BOOST_CHECK_EQUAL((*m2)(3, 0), 45.12);
251 BOOST_CHECK_EQUAL((*m2)(3, 1), 76.12);
252 BOOST_CHECK_EQUAL((*m2)(3, 2), 45.23);
253
254#ifdef NEKTAR_FULLDEBUG
255 BOOST_CHECK_THROW((*m2)(0, 3), ErrorUtil::NekError);
256 BOOST_CHECK_THROW((*m2)(4, 1), ErrorUtil::NekError);
257 BOOST_CHECK_THROW((*m2)(4, 3), ErrorUtil::NekError);
258#endif // NEKTAR_FULLDEBUG
259}
260
261BOOST_AUTO_TEST_CASE(TestDiagonalMatrixGetValue)
262{
263 UnitTests::RedirectCerrIfNeeded();
264 double data[] = {8.9, 3.4, 5.7};
265 Nektar::NekMatrix<double> m1(3, 3, data, eDIAGONAL);
266 std::shared_ptr<Nektar::Matrix<double>> m2(
267 new Nektar::NekMatrix<double>(3, 3, data, eDIAGONAL));
268 std::shared_ptr<Nektar::NekMatrix<double>> m3(
269 new Nektar::NekMatrix<double>(3, 3, data, eDIAGONAL));
270
271 BOOST_CHECK_EQUAL(m1(0, 0), 8.9);
272 BOOST_CHECK_EQUAL((*m2)(0, 0), 8.9);
273 BOOST_CHECK_EQUAL((*m3)(0, 0), 8.9);
274
275 BOOST_CHECK_EQUAL(m1(0, 1), 0.0);
276 BOOST_CHECK_EQUAL((*m2)(0, 1), 0.0);
277 BOOST_CHECK_EQUAL((*m3)(0, 1), 0.0);
278
279 BOOST_CHECK_EQUAL(m1(0, 2), 0.0);
280 BOOST_CHECK_EQUAL((*m2)(0, 2), 0.0);
281 BOOST_CHECK_EQUAL((*m3)(0, 2), 0.0);
282
283 BOOST_CHECK_EQUAL(m1(1, 0), 0.0);
284 BOOST_CHECK_EQUAL((*m2)(1, 0), 0.0);
285 BOOST_CHECK_EQUAL((*m3)(1, 0), 0.0);
286
287 BOOST_CHECK_EQUAL(m1(1, 1), 3.4);
288 BOOST_CHECK_EQUAL((*m2)(1, 1), 3.4);
289 BOOST_CHECK_EQUAL((*m3)(1, 1), 3.4);
290
291 BOOST_CHECK_EQUAL(m1(1, 2), 0.0);
292 BOOST_CHECK_EQUAL((*m2)(1, 2), 0.0);
293 BOOST_CHECK_EQUAL((*m3)(1, 2), 0.0);
294
295 BOOST_CHECK_EQUAL(m1(2, 0), 0.0);
296 BOOST_CHECK_EQUAL((*m2)(2, 0), 0.0);
297 BOOST_CHECK_EQUAL((*m3)(2, 0), 0.0);
298
299 BOOST_CHECK_EQUAL(m1(2, 1), 0.0);
300 BOOST_CHECK_EQUAL((*m2)(2, 1), 0.0);
301 BOOST_CHECK_EQUAL((*m3)(2, 1), 0.0);
302
303 BOOST_CHECK_EQUAL(m1(2, 2), 5.7);
304 BOOST_CHECK_EQUAL((*m2)(2, 2), 5.7);
305 BOOST_CHECK_EQUAL((*m3)(2, 2), 5.7);
306
307#ifdef NEKTAR_FULLDEBUG
308 BOOST_CHECK_THROW(m1(0, 4), ErrorUtil::NekError);
309 BOOST_CHECK_THROW((*m2)(0, 4), ErrorUtil::NekError);
310 BOOST_CHECK_THROW((*m3)(0, 4), ErrorUtil::NekError);
311
312 BOOST_CHECK_THROW(m1(4, 0), ErrorUtil::NekError);
313 BOOST_CHECK_THROW((*m2)(4, 0), ErrorUtil::NekError);
314 BOOST_CHECK_THROW((*m3)(4, 0), ErrorUtil::NekError);
315
316 BOOST_CHECK_THROW(m1(4, 4), ErrorUtil::NekError);
317 BOOST_CHECK_THROW((*m2)(4, 4), ErrorUtil::NekError);
318 BOOST_CHECK_THROW((*m3)(4, 4), ErrorUtil::NekError);
319#endif // NEKTAR_FULLDEBUG
320}
321
322BOOST_AUTO_TEST_CASE(TestFullNekMatrixSetValue)
323{
324 double data[] = {1.0, 3.0, 5.0, 2.0, 4.0, 6.0};
325 Nektar::NekMatrix<double> m1(3, 2, data);
326 std::shared_ptr<Nektar::Matrix<double>> m2(
327 new Nektar::NekMatrix<double>(3, 2, data));
328 std::shared_ptr<Nektar::NekMatrix<double>> m3(
329 new Nektar::NekMatrix<double>(3, 2, data));
330
331 m1.SetValue(0, 0, -1.0);
332 m2->SetValue(1, 1, -2.0);
333 m3->SetValue(2, 1, -3.0);
334
335 BOOST_CHECK_EQUAL(m1(0, 0), -1.0);
336 BOOST_CHECK_EQUAL(m1(0, 1), 2.0);
337 BOOST_CHECK_EQUAL(m1(1, 0), 3.0);
338 BOOST_CHECK_EQUAL(m1(1, 1), 4.0);
339 BOOST_CHECK_EQUAL(m1(2, 0), 5.0);
340 BOOST_CHECK_EQUAL(m1(2, 1), 6.0);
341
342 BOOST_CHECK_EQUAL((*m2)(0, 0), 1.0);
343 BOOST_CHECK_EQUAL((*m2)(0, 1), 2.0);
344 BOOST_CHECK_EQUAL((*m2)(1, 0), 3.0);
345 BOOST_CHECK_EQUAL((*m2)(1, 1), -2.0);
346 BOOST_CHECK_EQUAL((*m2)(2, 0), 5.0);
347 BOOST_CHECK_EQUAL((*m2)(2, 1), 6.0);
348
349 BOOST_CHECK_EQUAL((*m3)(0, 0), 1.0);
350 BOOST_CHECK_EQUAL((*m3)(0, 1), 2.0);
351 BOOST_CHECK_EQUAL((*m3)(1, 0), 3.0);
352 BOOST_CHECK_EQUAL((*m3)(1, 1), 4.0);
353 BOOST_CHECK_EQUAL((*m3)(2, 0), 5.0);
354 BOOST_CHECK_EQUAL((*m3)(2, 1), -3.0);
355}
356
357BOOST_AUTO_TEST_CASE(TestDiagonalNekMatrixSetValue)
358{
359 UnitTests::RedirectCerrIfNeeded();
360 double data[] = {8.9, 3.4, 5.7};
361 Nektar::NekMatrix<double> m1(3, 3, data, eDIAGONAL);
362 std::shared_ptr<Nektar::Matrix<double>> m2(
363 new Nektar::NekMatrix<double>(3, 3, data, eDIAGONAL));
364 std::shared_ptr<Nektar::NekMatrix<double>> m3(
365 new Nektar::NekMatrix<double>(3, 3, data, eDIAGONAL));
366
367 m1.SetValue(0, 0, 1.0);
368 m2->SetValue(1, 1, 2.0);
369 m3->SetValue(2, 2, 3.0);
370
371 BOOST_CHECK_EQUAL(m1.GetStorageType(), eDIAGONAL);
372 BOOST_CHECK_EQUAL(m2->GetStorageType(), eDIAGONAL);
373 BOOST_CHECK_EQUAL(m3->GetStorageType(), eDIAGONAL);
374
375 BOOST_CHECK_EQUAL(m1(0, 0), 1.0);
376 BOOST_CHECK_EQUAL(m1(1, 1), 3.4);
377 BOOST_CHECK_EQUAL(m1(2, 2), 5.7);
378
379 BOOST_CHECK_EQUAL((*m2)(0, 0), 8.9);
380 BOOST_CHECK_EQUAL((*m2)(1, 1), 2.0);
381 BOOST_CHECK_EQUAL((*m2)(2, 2), 5.7);
382
383 BOOST_CHECK_EQUAL((*m3)(0, 0), 8.9);
384 BOOST_CHECK_EQUAL((*m3)(1, 1), 3.4);
385 BOOST_CHECK_EQUAL((*m3)(2, 2), 3.0);
386
387#ifdef NEKTAR_FULLDEBUG
388 BOOST_CHECK_THROW(m1.SetValue(0, 3, 2.0), ErrorUtil::NekError);
389 BOOST_CHECK_THROW(m1.SetValue(3, 0, 2.0), ErrorUtil::NekError);
390 BOOST_CHECK_THROW(m1.SetValue(3, 3, 2.0), ErrorUtil::NekError);
391
392 BOOST_CHECK_THROW(m2->SetValue(0, 3, 2.0), ErrorUtil::NekError);
393 BOOST_CHECK_THROW(m2->SetValue(3, 0, 2.0), ErrorUtil::NekError);
394 BOOST_CHECK_THROW(m2->SetValue(3, 3, 2.0), ErrorUtil::NekError);
395
396 BOOST_CHECK_THROW(m3->SetValue(0, 3, 2.0), ErrorUtil::NekError);
397 BOOST_CHECK_THROW(m3->SetValue(3, 0, 2.0), ErrorUtil::NekError);
398 BOOST_CHECK_THROW(m3->SetValue(3, 3, 2.0), ErrorUtil::NekError);
399#endif // NEKTAR_FULLDEBUG
400}
401
402BOOST_AUTO_TEST_CASE(TestFullFullMatrixAddition)
403{
404 double lhs_data[] = {1.0, 4.0, 2.0, 5.0, 3.0, 6.0};
405 double rhs_data[] = {10.0, 13.0, 11.0, 14.0, 12.0, 15.0};
406
407 Nektar::NekMatrix<double> lhs(2, 3, lhs_data);
408 Nektar::NekMatrix<double> rhs(2, 3, rhs_data);
409 Nektar::NekMatrix<double> result = lhs + rhs;
410
411 BOOST_CHECK_EQUAL(result.GetRows(), 2);
412 BOOST_CHECK_EQUAL(result.GetColumns(), 3);
413 BOOST_CHECK_EQUAL(result.GetStorageType(), eFULL);
414 BOOST_CHECK_EQUAL(result(0, 0), 11.0);
415 BOOST_CHECK_EQUAL(result(0, 1), 13.0);
416 BOOST_CHECK_EQUAL(result(0, 2), 15.0);
417 BOOST_CHECK_EQUAL(result(1, 0), 17.0);
418 BOOST_CHECK_EQUAL(result(1, 1), 19.0);
419 BOOST_CHECK_EQUAL(result(1, 2), 21.0);
420}
421} // namespace MatrixUnitTests
422
423namespace UnitTests
424{
425
426using namespace Nektar;
427
428BOOST_AUTO_TEST_CASE(testNekMatrixConstruction)
429{
430 // Basic, dense matrix construction.
431 {
432 double buf[] = {1.0, 4.0, 7.0, 10.0, 2.0, 5.0,
433 8.0, 11.0, 3.0, 6.0, 9.0, 12.0};
434 NekMatrix<double> dynamic_matrix(4, 3, buf);
435
436 BOOST_CHECK(dynamic_matrix.GetRows() == 4);
437 BOOST_CHECK(dynamic_matrix.GetColumns() == 3);
438
439 for (unsigned int i = 0; i < 4; ++i)
440 {
441 for (unsigned int j = 0; j < 3; ++j)
442 {
443 BOOST_CHECK(dynamic_matrix(i, j) == buf[4 * j + i]);
444 }
445 }
446 }
447}
448
449BOOST_AUTO_TEST_CASE(testNekMatrixBasicMath)
450{
451 // Addition tests.
452 {
453 double buf[] = {1.0, 4.0, 7.0, 2.0, 5.0, 8.0, 3.0, 6.0, 9.0};
454 NekMatrix<double> m1(3, 3, buf);
455 NekMatrix<double> m2(3, 3, buf);
456 NekMatrix<double> m3 = m1 + m2;
457
458 for (unsigned int i = 0; i < 3; ++i)
459 {
460 for (unsigned int j = 0; j < 3; ++j)
461 {
462 BOOST_CHECK(m3(i, j) == buf[3 * j + i] + buf[3 * j + i]);
463 }
464 }
465
466 NekMatrix<double> m4(3, 3, buf);
467 NekMatrix<double> m5(3, 3, buf);
468 NekMatrix<double> m6 = m4 + m5;
469
470 for (unsigned int i = 0; i < 3; ++i)
471 {
472 for (unsigned int j = 0; j < 3; ++j)
473 {
474 BOOST_CHECK(m6(i, j) == buf[3 * j + i] + buf[3 * j + i]);
475 }
476 }
477 }
478 // Multiply
479
480 {
481 double buf1[] = {1, 4, 7, 2, 5, 8, 3, 6, 9};
482 double buf2[] = {10, 15, 19, 11, 16, 20, 12, 17, 21};
483 NekMatrix<double> lhs(3, 3, buf1);
484 NekMatrix<double> rhs(3, 3, buf2);
485
486 NekMatrix<double> result = lhs * rhs;
487
488 BOOST_CHECK(result.GetRows() == 3);
489 BOOST_CHECK(result.GetColumns() == 3);
490
491 double epsilon = 1e-12;
492 BOOST_CHECK_CLOSE(*result(0, 0), 97.0, epsilon);
493 BOOST_CHECK_CLOSE(*result(0, 1), 103.0, epsilon);
494 BOOST_CHECK_CLOSE(*result(0, 2), 109.0, epsilon);
495
496 BOOST_CHECK_CLOSE(*result(1, 0), 229.0, epsilon);
497 BOOST_CHECK_CLOSE(*result(1, 1), 244.0, epsilon);
498 BOOST_CHECK_CLOSE(*result(1, 2), 259.0, epsilon);
499
500 BOOST_CHECK_CLOSE(*result(2, 0), 361.0, epsilon);
501 BOOST_CHECK_CLOSE(*result(2, 1), 385.0, epsilon);
502 BOOST_CHECK_CLOSE(*result(2, 2), 409.0, epsilon);
503 }
504
505 {
506 double buf1[] = {1, 4, 7, 2, 5, 8, 3, 6, 9};
507 double buf2[] = {10, 15, 19, 11, 16, 20, 12, 17, 21, 14, 18, 22};
508 NekMatrix<double> lhs(3, 3, buf1);
509 NekMatrix<double> rhs(3, 4, buf2);
510
511 NekMatrix<double> result = lhs * rhs;
512
513 BOOST_CHECK(result.GetRows() == 3);
514 BOOST_CHECK(result.GetColumns() == 4);
515
516 double epsilon = 1e-12;
517 BOOST_CHECK_CLOSE(*result(0, 0), 97.0, epsilon);
518 BOOST_CHECK_CLOSE(*result(0, 1), 103.0, epsilon);
519 BOOST_CHECK_CLOSE(*result(0, 2), 109.0, epsilon);
520 BOOST_CHECK_CLOSE(*result(0, 3), 116.0, epsilon);
521
522 BOOST_CHECK_CLOSE(*result(1, 0), 229.0, epsilon);
523 BOOST_CHECK_CLOSE(*result(1, 1), 244.0, epsilon);
524 BOOST_CHECK_CLOSE(*result(1, 2), 259.0, epsilon);
525 BOOST_CHECK_CLOSE(*result(1, 3), 278.0, epsilon);
526
527 BOOST_CHECK_CLOSE(*result(2, 0), 361.0, epsilon);
528 BOOST_CHECK_CLOSE(*result(2, 1), 385.0, epsilon);
529 BOOST_CHECK_CLOSE(*result(2, 2), 409.0, epsilon);
530 BOOST_CHECK_CLOSE(*result(2, 3), 440.0, epsilon);
531 }
532
533 //
534 // // Transpose
535 //
536 // // Determinant.
537 //
538 // // Invert/Check for singularity.
539 //
540 // // Eigenvalues/vectors
541 //
542 // // Condition number wrt various norms.
543 //
544 // // Various norm computations.
545 //
546 // // LU Decomposition? More appropriate in LinAlg?
547}
548
549} // namespace UnitTests
550} // namespace Nektar
Nektar::MatrixUnitTests::BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(TestFullMatrixInversion)
Definition testNekMatrix.cpp:45
Nektar::UnitTests::RedirectCerrIfNeeded
void RedirectCerrIfNeeded()
Definition util.cpp:41
Nektar::BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(TestCanGetRawPtr)
Definition TestCanGetRawPtr.cpp:44
Generated on Sat Jul 4 2026 01:02:24 for Nektar++ by doxygen 1.9.8