Nektar++
TestBandedMatrixOperations.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: TestBandedMatrixOperations.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 <LibUtilities/LinearAlgebra/NekLinSys.hpp>
36#include <LibUtilities/LinearAlgebra/NekMatrix.hpp>
37#include <LibUtilities/LinearAlgebra/NekTypeDefs.hpp>
38#include <UnitTests/LibUtilities/LinearAlgebra/TestCombinationRunner.h>
39
40#include <boost/test/tools/floating_point_comparison.hpp>
41#include <boost/test/unit_test.hpp>
42
43namespace Nektar
44{
45namespace BandedMatrixVectorMultiplicationUnitTests
46{
47BOOST_AUTO_TEST_CASE(TestDirectBlasCall)
48{
49 // [1 2 0 0 ]
50 // [3 4 5 0 ]
51 // [0 6 7 8 ]
52 // [0 0 9 10]
53 NekDouble x[] = {1, 2, 3, 4};
54 NekDouble y[] = {0, 0, 0, 0};
55
56 NekDouble a[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0};
57 Blas::Dgbmv('T', 4, 4, 1, 1, 1.0, a, 3, x, 1, 0.0, y, 1);
58
59 NekDouble expected_result_buf[] = {5, 26, 65, 67};
60 NekVector<NekDouble> expected_result(4, expected_result_buf);
61 NekVector<NekDouble> result(4, y);
62 BOOST_CHECK_EQUAL(expected_result, result);
63}
64
65BOOST_AUTO_TEST_CASE(TestSquareDirectBlasCall)
66{
67 // [1 2 3 4 ]
68 // [5 6 7 8 ]
69 // [9 10 11 12]
70 // [13 14 15 16]
71 NekDouble x[] = {1, 2, 3, 4};
72 NekDouble y[] = {0, 0, 0, 0};
73
74 NekDouble a[] = {0, 0, 0, 1, 2, 3, 4, 0, 0, 5, 6, 7, 8, 0,
75 0, 9, 10, 11, 12, 0, 0, 13, 14, 15, 16, 0, 0, 0};
76 Blas::Dgbmv('T', 4, 4, 3, 3, 1.0, a, 7, x, 1, 0.0, y, 1);
77
78 NekDouble expected_result_buf[] = {30, 70, 110, 150};
79 NekVector<NekDouble> expected_result(4, expected_result_buf);
80 NekVector<NekDouble> result(4, y);
81 BOOST_CHECK_EQUAL(expected_result, result);
82}
83
84BOOST_AUTO_TEST_CASE(TestDifferentNumberOfLowerAndUpperDiagonalsDirectBlasCall)
85{
86 // [1 2 11 0 ]
87 // [3 4 5 12 ]
88 // [0 6 7 8 ]
89 // [0 0 9 10]
90 NekDouble x[] = {1, 2, 3, 4};
91 NekDouble y[] = {0, 0, 0, 0};
92
93 NekDouble a[] = {0, 1, 2, 11, 3, 4, 5, 12, 6, 7, 8, 0, 9, 10, 0, 0};
94
95 Blas::Dgbmv('T', 4, 4, 2, 1, 1.0, a, 4, x, 1, 0.0, y, 1);
96
97 NekDouble expected_result_buf[] = {38, 74, 65, 67};
98 NekVector<NekDouble> expected_result(4, expected_result_buf);
99 NekVector<NekDouble> result(4, y);
100 BOOST_CHECK_EQUAL(expected_result, result);
101}
102
103BOOST_AUTO_TEST_CASE(TestLargerPackedMatrixThanNormal)
104{
105 // [ 1 2 0 ]
106 // [ 3 4 5 ]
107 // [ 6 7 8 ]
108 NekDouble x[] = {1, 2, 3};
109 NekDouble y[] = {0, 0, 0};
110
111 NekDouble a[] = {0, 1, 3, 6, 2, 4, 7, 0, 5, 8, 0, 0};
112
113 Blas::Dgbmv('N', 3, 3, 2, 1, 1.0, a, 4, x, 1, 0.0, y, 1);
114
115 NekDouble expected_result_buf[] = {5, 26, 44};
116 NekVector<NekDouble> expected_result(3, expected_result_buf);
117 NekVector<NekDouble> result(3, y);
118 BOOST_CHECK_EQUAL(expected_result, result);
119}
120
121BOOST_AUTO_TEST_CASE(TestStandardMatrixVectorMultiply)
122{
123
124 {
125 // This is the matrix
126 // [ 1 2 0 0 ]
127 // [ 3 4 5 0 ]
128 // [ 0 6 7 8 ]
129 // [ 0 0 9 10]
130
131 NekDouble buf[] = {0, 1, 3, 2, 4, 6, 5, 7, 9, 8, 10, 0};
132
133 std::shared_ptr<DenseMatrix> m(
134 new DenseMatrix(4, 4, buf, eBANDED, 1, 1));
135
136 std::shared_ptr<ScaledMatrix> scaled(new ScaledMatrix(2.0, m));
137
138 NekDouble vector_buf[] = {1.0, 2.0, 3.0, 4.0};
139 NekVector<NekDouble> v(4, vector_buf);
140
141 NekVector<NekDouble> result = (*m) * v;
142 NekVector<NekDouble> scaledResult = (*scaled) * v;
143
144 NekDouble expected_result_buf[] = {5, 26, 65, 67};
145 NekVector<NekDouble> expected_result(4, expected_result_buf);
146 NekVector<NekDouble> expectedScaledResult = 2.0 * expected_result;
147 BOOST_CHECK_EQUAL(expected_result, result);
148
149 BOOST_CHECK_EQUAL(expectedScaledResult, scaledResult);
150 }
151}
152
153BOOST_AUTO_TEST_CASE(
154 TestUnequalNumbersOfSubAndSuperDiagonalsMatrixVectorMultiply)
155{
156 {
157 typedef NekMatrix<NekDouble> MatrixType;
158
159 // This is the matrix
160 // [ 1 6 10 0 0 ]
161 // [ 13 2 7 11 0 ]
162 // [ 0 14 3 8 12 ]
163 // [ 0 0 15 4 9 ]
164 // [ 0 0 0 16 5 ]
165
166 NekDouble buf[] = {0, 0, 1, 13, 0, 6, 2, 14, 10, 7,
167 3, 15, 11, 8, 4, 16, 12, 9, 5, 0};
168
169 MatrixType m(5, 5, buf, eBANDED, 1, 2);
170
171 NekDouble vector_buf[] = {1, 2, 3, 4, 5};
172 NekVector<NekDouble> v(5, vector_buf);
173
174 NekVector<NekDouble> result = m * v;
175
176 NekDouble expected_result_buf[] = {43, 82, 129, 106, 89};
177 NekVector<NekDouble> expected_result(5, expected_result_buf);
178 BOOST_CHECK_EQUAL(expected_result, result);
179 }
180}
181
182BOOST_AUTO_TEST_CASE(TestNoSubdiagonalsTwoSuperDiagonalsMatrixVectorMultiply)
183{
184 {
185 typedef NekMatrix<NekDouble> MatrixType;
186
187 // This is the matrix
188 // [ 1 6 10 0 0 ]
189 // [ 0 2 7 11 0 ]
190 // [ 0 0 3 8 12 ]
191 // [ 0 0 0 4 9 ]
192 // [ 0 0 0 0 5 ]
193
194 NekDouble buf[] = {0, 0, 1, 0, 6, 2, 10, 7, 3, 11, 8, 4, 12, 9, 5};
195
196 MatrixType m(5, 5, buf, eBANDED, 0, 2);
197
198 NekDouble vector_buf[] = {1, 2, 3, 4, 5};
199 NekVector<NekDouble> v(5, vector_buf);
200
201 NekVector<NekDouble> result = m * v;
202
203 NekDouble expected_result_buf[] = {43, 69, 101, 61, 25};
204 NekVector<NekDouble> expected_result(5, expected_result_buf);
205 BOOST_CHECK_EQUAL(expected_result, result);
206 }
207}
208
209BOOST_AUTO_TEST_CASE(TestNoSubdiagonalsThreeSuperDiagonalsMatrixVectorMultiply)
210{
211 {
212 typedef NekMatrix<NekDouble> MatrixType;
213
214 // This is the matrix
215 // [ 1 6 10 2 0 ]
216 // [ 0 2 7 11 5 ]
217 // [ 0 0 3 8 12 ]
218 // [ 0 0 0 4 9 ]
219 // [ 0 0 0 0 5 ]
220
221 NekDouble buf[] = {0, 0, 0, 1, 0, 0, 6, 2, 0, 10,
222 7, 3, 2, 11, 8, 4, 5, 12, 9, 5};
223
224 MatrixType m(5, 5, buf, eBANDED, 0, 3);
225
226 NekDouble vector_buf[] = {1, 2, 3, 4, 5};
227 NekVector<NekDouble> v(5, vector_buf);
228
229 NekVector<NekDouble> result = m * v;
230
231 NekDouble expected_result_buf[] = {51, 94, 101, 61, 25};
232 NekVector<NekDouble> expected_result(5, expected_result_buf);
233 BOOST_CHECK_EQUAL(expected_result, result);
234 }
235}
236
237BOOST_AUTO_TEST_CASE(TestNoSubdiagonalsFourSuperDiagonalsMatrixVectorMultiply)
238{
239 {
240 typedef NekMatrix<NekDouble> MatrixType;
241
242 // This is the matrix
243 // [ 1 6 10 2 8 ]
244 // [ 0 2 7 11 5 ]
245 // [ 0 0 3 8 12 ]
246 // [ 0 0 0 4 9 ]
247 // [ 0 0 0 0 5 ]
248
249 NekDouble buf[] = {0, 0, 0, 0, 1, 0, 0, 0, 6, 2, 0, 0, 10,
250 7, 3, 0, 2, 11, 8, 4, 8, 5, 12, 9, 5};
251
252 MatrixType m(5, 5, buf, eBANDED, 0, 4);
253
254 NekDouble vector_buf[] = {1, 2, 3, 4, 5};
255 NekVector<NekDouble> v(5, vector_buf);
256
257 NekVector<NekDouble> result = m * v;
258
259 NekDouble expected_result_buf[] = {91, 94, 101, 61, 25};
260 NekVector<NekDouble> expected_result(5, expected_result_buf);
261 BOOST_CHECK_EQUAL(expected_result, result);
262 }
263}
264
265BOOST_AUTO_TEST_CASE(Test2Sub2SuperMatrixVectorMultiply)
266{
267 {
268 typedef NekMatrix<NekDouble> MatrixType;
269
270 // This is the matrix
271 // [ 1 6 10 0 0 ]
272 // [ 13 2 7 11 0 ]
273 // [ 7 14 3 8 12 ]
274 // [ 0 5 15 4 9 ]
275 // [ 0 0 4 16 5 ]
276
277 NekDouble buf[] = {0, 0, 1, 13, 7, 0, 6, 2, 14, 5, 10, 7, 3,
278 15, 4, 11, 8, 4, 16, 0, 12, 9, 5, 0, 0};
279
280 MatrixType m(5, 5, buf, eBANDED, 2, 2);
281
282 NekDouble vector_buf[] = {1, 2, 3, 4, 5};
283 NekVector<NekDouble> v(5, vector_buf);
284
285 NekVector<NekDouble> result = m * v;
286
287 NekDouble expected_result_buf[] = {43, 82, 136, 116, 101};
288 NekVector<NekDouble> expected_result(5, expected_result_buf);
289 BOOST_CHECK_EQUAL(expected_result, result);
290 }
291}
292
293BOOST_AUTO_TEST_CASE(Test3Sub2SuperMatrixVectorMultiply)
294{
295 {
296 typedef NekMatrix<NekDouble> MatrixType;
297
298 // This is the matrix
299 // [ 1 6 10 0 0 ]
300 // [ 13 2 7 11 0 ]
301 // [ 7 14 3 8 12 ]
302 // [ 20 5 15 4 9 ]
303 // [ 0 21 4 16 5 ]
304
305 NekDouble buf[] = {0, 0, 1, 13, 7, 20, 0, 6, 2, 14, 5, 21, 10, 7, 3,
306 15, 4, 0, 11, 8, 4, 16, 0, 0, 12, 9, 5, 0, 0, 0};
307
308 MatrixType m(5, 5, buf, eBANDED, 3, 2);
309
310 NekDouble vector_buf[] = {1, 2, 3, 4, 5};
311 NekVector<NekDouble> v(5, vector_buf);
312
313 NekVector<NekDouble> result = m * v;
314
315 NekDouble expected_result_buf[] = {43, 82, 136, 136, 143};
316 NekVector<NekDouble> expected_result(5, expected_result_buf);
317 BOOST_CHECK_EQUAL(expected_result, result);
318 }
319}
320
321BOOST_AUTO_TEST_CASE(Test4Sub2SuperMatrixVectorMultiply)
322{
323 {
324 typedef NekMatrix<NekDouble> MatrixType;
325
326 // This is the matrix
327 // [ 1 6 10 0 0 ]
328 // [ 13 2 7 11 0 ]
329 // [ 7 14 3 8 12 ]
330 // [ 20 5 15 4 9 ]
331 // [ 30 21 4 16 5 ]
332
333 NekDouble buf[] = {0, 0, 1, 13, 7, 20, 30, 0, 6, 2, 14, 5,
334 21, 0, 10, 7, 3, 15, 4, 0, 0, 11, 8, 4,
335 16, 0, 0, 0, 12, 9, 5, 0, 0, 0, 0};
336
337 MatrixType m(5, 5, buf, eBANDED, 4, 2);
338
339 NekDouble vector_buf[] = {1, 2, 3, 4, 5};
340 NekVector<NekDouble> v(5, vector_buf);
341
342 NekVector<NekDouble> result = m * v;
343
344 NekDouble expected_result_buf[] = {43, 82, 136, 136, 173};
345 NekVector<NekDouble> expected_result(5, expected_result_buf);
346 BOOST_CHECK_EQUAL(expected_result, result);
347 }
348}
349
350BOOST_AUTO_TEST_CASE(TestDiagonalOnlyMatrixVectorMultiply)
351{
352 {
353 typedef NekMatrix<NekDouble> MatrixType;
354
355 // This is the matrix
356 // [ 5 0 0 0 0 ]
357 // [ 0 1 0 0 0 ]
358 // [ 0 0 9 0 0 ]
359 // [ 0 0 0 4 0 ]
360 // [ 0 0 0 0 2 ]
361
362 NekDouble buf[] = {5, 1, 9, 4, 2};
363
364 MatrixType m(5, 5, buf, eBANDED, 0, 0);
365
366 NekDouble vector_buf[] = {1, 2, 3, 4, 5};
367 NekVector<NekDouble> v(5, vector_buf);
368
369 NekVector<NekDouble> result = m * v;
370
371 NekDouble expected_result_buf[] = {5, 2, 27, 16, 10};
372 NekVector<NekDouble> expected_result(5, expected_result_buf);
373 BOOST_CHECK_EQUAL(expected_result, result);
374 }
375}
376
377BOOST_AUTO_TEST_CASE(TestBandedMatrixLinearSystemSolves)
378{
379 typedef NekMatrix<double> MatrixType;
380
381 {
382 double buf[] = {0, 0, -.23, -6.98, 0, 2.54, 2.46, 2.56,
383 -3.66, -2.73, 2.46, -4.78, -2.13, 4.07, -3.82, 0};
384
385 MatrixType m(4, 4, buf, eBANDED, 1, 2);
386 LinearSystem l(m);
387
388 double b_buf[] = {4.42, 27.13, -6.14, 10.5};
389 NekVector<double> b(4, b_buf);
390
391 NekVector<double> result = l.Solve(b);
392 }
393
394 {
395
396 // This is the matrix
397 // [ 1 6 10 0 0 ]
398 // [ 13 2 7 11 0 ]
399 // [ 7 14 3 8 12 ]
400 // [ 20 5 15 4 9 ]
401 // [ 30 21 4 16 5 ]
402
403 double buf[] = {0, 0, 1, 13, 7, 20, 30, 0, 6, 2, 14, 5,
404 21, 0, 10, 7, 3, 15, 4, 0, 0, 11, 8, 4,
405 16, 0, 0, 0, 12, 9, 5, 0, 0, 0, 0};
406
407 MatrixType m(5, 5, buf, eBANDED, 4, 2);
408
409 double b_buf[] = {43, 82, 136, 136, 173};
410 NekVector<double> b(5, b_buf);
411
412 LinearSystem l(m);
413 NekVector<double> result(5, 0.0);
414 result = l.Solve(b);
415
416 double expected_result_buf[] = {1, 2, 3, 4, 5};
417 NekVector<double> expected_result(5, expected_result_buf);
418
419 BOOST_CHECK_CLOSE(expected_result[0], result[0], .00001);
420 BOOST_CHECK_CLOSE(expected_result[1], result[1], .00001);
421 BOOST_CHECK_CLOSE(expected_result[2], result[2], .00001);
422 BOOST_CHECK_CLOSE(expected_result[3], result[3], .00001);
423 BOOST_CHECK_CLOSE(expected_result[4], result[4], .00001);
424 }
425}
426
427} // namespace BandedMatrixVectorMultiplicationUnitTests
428
429namespace BandedMatrixMatrixMultiplicationTests
430{
431BOOST_AUTO_TEST_CASE(TestMatrixMatrixMultiplication)
432{
433 // [ 2 10 0 0 ] [ 30 38 0 0 ]
434 // [16 4 12 0 ] [ 0 32 40 0 ]
435 // [22 18 6 14 ] * [ 0 0 34 42 ]
436 // [ 0 24 20 8 ] [ 0 0 0 36 ]
437
438 {
439 double lhs_buf[] = {0, 2, 16, 22, 10, 4, 18, 24,
440 12, 6, 20, 0, 14, 8, 0, 0};
441 double rhs_buf[] = {0, 30, 38, 32, 40, 34, 42, 36};
442
443 NekMatrix<double> lhs1(4, 4, lhs_buf, eBANDED, 2, 1);
444 std::shared_ptr<
445 NekMatrix<NekMatrix<double, StandardMatrixTag>, ScaledMatrixTag>>
446 lhs2;
447 std::shared_ptr<NekMatrix<NekMatrix<double>, BlockMatrixTag>> lhs3;
448
449 NekMatrix<double> rhs1(4, 4, rhs_buf, eBANDED, 0, 1);
450 std::shared_ptr<NekMatrix<NekMatrix<double>, ScaledMatrixTag>> rhs2;
451 std::shared_ptr<NekMatrix<NekMatrix<double>, BlockMatrixTag>> rhs3;
452
453 GenerateMatrices(lhs1, 2.0, 2, 2, lhs2, lhs3);
454 GenerateMatrices(rhs1, 2.0, 2, 2, rhs2, rhs3);
455
456 double result_buf[] = {60, 480, 660, 0, 396, 736, 1412, 768,
457 400, 568, 924, 1640, 0, 504, 756, 1128};
458 NekMatrix<double> result(4, 4, result_buf);
459 }
460}
461} // namespace BandedMatrixMatrixMultiplicationTests
462} // namespace Nektar
Nektar::LinearSystem::Solve
RawType_t< VectorType > Solve(const VectorType &b)
Definition: NekLinSys.hpp:448
Blas::Dgbmv
static void Dgbmv(const char &trans, const int &m, const int &n, const int &kl, const int &ku, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = alpha A x plus beta y where A[m x n] is banded.
Definition: Blas.hpp:243
Nektar::GenerateMatrices
void GenerateMatrices(const NekMatrix< NumberType, StandardMatrixTag > &m1, NumberType scale, unsigned int blockRows, unsigned int blockColumns, std::shared_ptr< NekMatrix< NekMatrix< NumberType, StandardMatrixTag >, ScaledMatrixTag > > &m2, std::shared_ptr< NekMatrix< NekMatrix< NumberType >, BlockMatrixTag > > &m3)
Definition: TestCombinationRunner.h:98
Nektar::ScaledMatrix
NekMatrix< DenseMatrix, ScaledMatrixTag > ScaledMatrix
Definition: NekTypeDefs.hpp:55
Nektar::DenseMatrix
NekMatrix< NekDouble, StandardMatrixTag > DenseMatrix
Definition: NekTypeDefs.hpp:51
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Generated on Wed Jan 31 2024 09:35:57 for Nektar++ by doxygen 1.9.4