Nektar++
TestFullMatrixOperations.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: TestFullMatrixOperations.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 <UnitTests/CountedObject.h>
39
40#include <boost/test/tools/floating_point_comparison.hpp>
41#include <boost/test/unit_test.hpp>
42
43namespace Nektar::FullMatrixOperationsUnitTests
44{
45template <typename DataType, typename MatrixType>
46int foo([[maybe_unused]] NekMatrix<DataType, MatrixType> &d)
47{
48 return 1;
49}
50
51template <typename DataType>
52int foo([[maybe_unused]] NekMatrix<DataType, BlockMatrixTag> &d)
53{
54 return 2;
55}
56
57BOOST_AUTO_TEST_CASE(TestDoubleSquareFullVectorMultiplication)
58{
59 NekMatrix<double> m1;
60 NekMatrix<NekMatrix<double>, BlockMatrixTag> m2(2, 2, 2, 2);
61 foo(m1);
62 foo(m2);
63 double m_buf[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
64 double v_buf[] = {4, 5, 6};
65
66 NekMatrix<double> m(3, 3, m_buf);
67 NekVector<double> v(3, v_buf);
68
69 double expected_result_buf[] = {66, 81, 96};
70 NekVector<double> expected_result(3, expected_result_buf);
71
72 NekVector<double> result = m * v;
73
74 BOOST_CHECK_EQUAL(expected_result, result);
75}
76
77BOOST_AUTO_TEST_CASE(TestDoubleSquareFullVectorMultiplicationWithAliasing)
78{
79 double m_buf[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
80 double v_buf[] = {4, 5, 6};
81
82 NekMatrix<double> m(3, 3, m_buf);
83 NekVector<double> v(3, v_buf);
84
85 double expected_result_buf[] = {66, 81, 96};
86 NekVector<double> expected_result(3, expected_result_buf);
87
88 v = m * v;
89
90 BOOST_CHECK_EQUAL(expected_result, v);
91}
92
93BOOST_AUTO_TEST_CASE(
94 TestDoubleSquareFullVectorMultiplicationWithSharedArrayAliasing)
95{
96 double m_buf[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
97 double v_buf[] = {4, 5, 6};
98
99 Array<OneD, double> vector_array(3, v_buf);
100
101 NekMatrix<double> m(3, 3, m_buf);
102 NekVector<double> v(3, vector_array, eWrapper);
103 NekVector<double> result(3, vector_array, eWrapper);
104
105 double expected_result_buf[] = {66, 81, 96};
106 NekVector<double> expected_result(3, expected_result_buf);
107
108 result = m * v;
109
110 BOOST_CHECK_EQUAL(expected_result, result);
111 BOOST_CHECK_EQUAL(vector_array[0], 66);
112 BOOST_CHECK_EQUAL(vector_array[1], 81);
113 BOOST_CHECK_EQUAL(vector_array[2], 96);
114}
115
116BOOST_AUTO_TEST_CASE(
117 TestDoubleSquareFullVectorMultiplicationWithSharedArrayAliasingAndMatrixTranspose)
118{
119 double m_buf[] = {1, 4, 7, 2, 5, 8, 3, 6, 9};
120 double v_buf[] = {4, 5, 6};
121
122 Array<OneD, double> vector_array(3, v_buf);
123
124 NekMatrix<double> m(3, 3, m_buf);
125 NekVector<double> v(3, vector_array, eWrapper);
126 NekVector<double> result(3, vector_array, eWrapper);
127
128 double expected_result_buf[] = {66, 81, 96};
129 NekVector<double> expected_result(3, expected_result_buf);
130
131 result = Transpose(m) * v;
132
133 BOOST_CHECK_EQUAL(expected_result, result);
134 BOOST_CHECK_EQUAL(vector_array[0], 66);
135 BOOST_CHECK_EQUAL(vector_array[1], 81);
136 BOOST_CHECK_EQUAL(vector_array[2], 96);
137}
138
139BOOST_AUTO_TEST_CASE(
140 TestDoubleSquareFullVectorMultiplicationWithSharedArrayAliasingAndOverlap)
141{
142 double m_buf[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
143 double v_buf[] = {4, 5, 6, 0};
144
145 Array<OneD, double> vector_array(4, v_buf);
146 Array<OneD, double> tmp;
147
148 NekMatrix<double> m(3, 3, m_buf);
149 NekVector<double> v(3, vector_array, eWrapper);
150 NekVector<double> result(3, tmp = vector_array + 1, eWrapper);
151
152 double expected_result_buf[] = {66, 81, 96};
153 NekVector<double> expected_result(3, expected_result_buf);
154
155 result = m * v;
156
157 BOOST_CHECK_EQUAL(expected_result, result);
158 BOOST_CHECK_EQUAL(vector_array[1], 66);
159 BOOST_CHECK_EQUAL(vector_array[2], 81);
160 BOOST_CHECK_EQUAL(vector_array[3], 96);
161}
162
163BOOST_AUTO_TEST_CASE(TestScaledDoubleSquareFullVectorMultiplication)
164{
165 double m_buf[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
166 double v_buf[] = {4, 5, 6};
167
168 std::shared_ptr<NekMatrix<double>> inner(
169 new NekMatrix<double>(3, 3, m_buf));
170 NekMatrix<NekMatrix<double>, ScaledMatrixTag> m(7.0, inner);
171 NekVector<double> v(3, v_buf);
172
173 double expected_result_buf[] = {462, 567, 672};
174 NekVector<double> expected_result(3, expected_result_buf);
175
176 NekVector<double> result = m * v;
177
178 BOOST_CHECK_EQUAL(expected_result, result);
179}
180
181BOOST_AUTO_TEST_CASE(TestThreeMatrixMultiplication)
182{
183 double buf1[] = {1.0, 2.0, 3.0, 4.0};
184 double buf2[] = {5.0, 6.0, 7.0, 8.0};
185 double buf3[] = {9.0, 10.0, 11.0, 12.0};
186
187 NekMatrix<double> m1(2, 2, buf1);
188 NekMatrix<double> m2(2, 2, buf2);
189 NekMatrix<double> m3(2, 2, buf3);
190
191 double expected_result_buf[] = {517, 766, 625, 926};
192 NekMatrix<double> expected_result(2, 2, expected_result_buf);
193
194 NekMatrix<double> result1 = m1 * m2 * m3;
195 NekMatrix<double> result2 = m1 * (m2 * m3);
196 BOOST_CHECK_EQUAL(expected_result, m1 * m2 * m3);
197
198 BOOST_CHECK_EQUAL(expected_result, m1 * (m2 * m3));
199}
200
201BOOST_AUTO_TEST_CASE(TestThreeMatrixMultiplicationWithTranspose)
202{
203 double buf1[] = {1.0, 2.0, 3.0, 4.0};
204 double buf2[] = {5.0, 6.0, 7.0, 8.0};
205 double buf3[] = {9.0, 11.0, 10.0, 12.0};
206
207 NekMatrix<double> m1(2, 2, buf1);
208 NekMatrix<double> m2(2, 2, buf2);
209 NekMatrix<double> m3(2, 2, buf3);
210
211 double expected_result_buf[] = {517, 766, 625, 926};
212 NekMatrix<double> expected_result(2, 2, expected_result_buf);
213
214 BOOST_CHECK_EQUAL(expected_result, m1 * m2 * Transpose(m3));
215}
216
217BOOST_AUTO_TEST_CASE(TestThreeWrappedMatrixMultiplicationWithTranspose)
218{
219 double buf1[] = {1.0, 2.0, 3.0, 4.0};
220 double buf2[] = {5.0, 6.0, 7.0, 8.0};
221 double buf3[] = {9.0, 11.0, 10.0, 12.0};
222 double out_buf[] = {0.0, 0.0, 0.0, 0.0};
223
224 Array<OneD, double> array_buf1(4, buf1);
225 Array<OneD, double> array_buf2(4, buf2);
226 Array<OneD, double> array_buf3(4, buf3);
227 Array<OneD, double> array_out_buf(4, out_buf);
228
229 NekMatrix<double> m1(2, 2, array_buf1, eWrapper);
230 NekMatrix<double> m2(2, 2, array_buf2, eWrapper);
231 NekMatrix<double> m3(2, 2, array_buf3, eWrapper);
232 NekMatrix<double> result(2, 2, array_out_buf, eWrapper);
233
234 double expected_result_buf[] = {517, 766, 625, 926};
235 NekMatrix<double> expected_result(2, 2, expected_result_buf);
236
237 result = m1 * m2 * Transpose(m3);
238 BOOST_CHECK_EQUAL(expected_result, result);
239 BOOST_CHECK_EQUAL(array_out_buf[0], 517.0);
240 BOOST_CHECK_EQUAL(array_out_buf[1], 766.0);
241 BOOST_CHECK_EQUAL(array_out_buf[2], 625.0);
242 BOOST_CHECK_EQUAL(array_out_buf[3], 926.0);
243 BOOST_CHECK_EQUAL(result(0, 0), 517.0);
244 BOOST_CHECK_EQUAL(result(1, 0), 766.0);
245 BOOST_CHECK_EQUAL(result(0, 1), 625.0);
246 BOOST_CHECK_EQUAL(result(1, 1), 926.0);
247}
248} // namespace Nektar::FullMatrixOperationsUnitTests
Nektar::FullMatrixOperationsUnitTests::BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(TestDoubleSquareFullVectorMultiplication)
Definition: TestFullMatrixOperations.cpp:57
Nektar::FullMatrixOperationsUnitTests::foo
int foo(NekMatrix< DataType, MatrixType > &d)
Definition: TestFullMatrixOperations.cpp:46
Nektar::UnitTests::d
std::vector< double > d(NPUPPER *NPUPPER)
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:235
Generated on Wed Jan 31 2024 09:35:57 for Nektar++ by doxygen 1.9.4