Nektar++
TestStandardMatrix.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: TestStandardMatrix.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#define BOOST_TEST_MODULE LinearAlgebraUnitTests test
36#include <boost/test/unit_test.hpp>
37
39#include <boost/test/tools/floating_point_comparison.hpp>
40#include <boost/test/unit_test.hpp>
41
43#include <functional>
44
45namespace Nektar
46{
47namespace StandardMatrixUnitTests
48{
50{
51 double buf1[] = {1, 4, 2, 5, 3, 6};
52
53 NekMatrix<double> m1(2, 3, buf1);
54 const NekMatrix<double> &m2 = m1;
55
56 NekMatrix<double>::iterator m1_iter = m1.begin();
57 NekMatrix<double>::const_iterator m2_iter = m2.begin();
58
59 BOOST_CHECK_EQUAL(1, *m1_iter);
60 BOOST_CHECK_EQUAL(1, *m2_iter);
61 BOOST_CHECK_EQUAL(1, *m1_iter++);
62 BOOST_CHECK_EQUAL(1, *m2_iter++);
63 BOOST_CHECK_EQUAL(4, *m1_iter++);
64 BOOST_CHECK_EQUAL(4, *m2_iter++);
65 BOOST_CHECK_EQUAL(2, *m1_iter++);
66 BOOST_CHECK_EQUAL(2, *m2_iter++);
67 BOOST_CHECK_EQUAL(5, *m1_iter++);
68 BOOST_CHECK_EQUAL(5, *m2_iter++);
69 BOOST_CHECK_EQUAL(3, *m1_iter++);
70 BOOST_CHECK_EQUAL(3, *m2_iter++);
71 BOOST_CHECK_EQUAL(6, *m1_iter++);
72 BOOST_CHECK_EQUAL(6, *m2_iter++);
73 BOOST_CHECK(m1_iter == m1.end());
74 BOOST_CHECK(m2_iter == m2.end());
75
76 m1.Transpose();
77 m1_iter = m1.begin();
78 m2_iter = m2.begin();
79
80 BOOST_CHECK_EQUAL(1, *m1_iter);
81 BOOST_CHECK_EQUAL(1, *m2_iter);
82 BOOST_CHECK_EQUAL(2, *(++m1_iter));
83 BOOST_CHECK_EQUAL(2, *(++m2_iter));
84 BOOST_CHECK_EQUAL(3, *(++m1_iter));
85 BOOST_CHECK_EQUAL(3, *(++m2_iter));
86 BOOST_CHECK_EQUAL(4, *(++m1_iter));
87 BOOST_CHECK_EQUAL(4, *(++m2_iter));
88 BOOST_CHECK_EQUAL(5, *(++m1_iter));
89 BOOST_CHECK_EQUAL(5, *(++m2_iter));
90 BOOST_CHECK_EQUAL(6, *(++m1_iter));
91 BOOST_CHECK_EQUAL(6, *(++m2_iter));
92 BOOST_CHECK(m1_iter != m1.end());
93 BOOST_CHECK(m2_iter != m2.end());
94 ++m1_iter;
95 ++m2_iter;
96 BOOST_CHECK(m1_iter == m1.end());
97 BOOST_CHECK(m2_iter == m2.end());
98}
99
100BOOST_AUTO_TEST_CASE(TestWrappedCopyConstructor)
101{
102 double buf1[] = {1, 4, 2, 5, 3, 6};
103
104 NekMatrix<double> m1(2, 3, buf1);
106
107 BOOST_CHECK_EQUAL(1.0, m1(0, 0));
108 BOOST_CHECK_EQUAL(1.0, m2(0, 0));
109
110 m1(0, 0) = -4.5;
111
112 BOOST_CHECK_EQUAL(-4.5, m1(0, 0));
113 BOOST_CHECK_EQUAL(-4.5, m2(0, 0));
114
115 BOOST_CHECK_EQUAL(4.0, m1(1, 0));
116 BOOST_CHECK_EQUAL(4.0, m2(1, 0));
117
118 m2(1, 0) = 6.7;
119
120 BOOST_CHECK_EQUAL(6.7, m1(1, 0));
121 BOOST_CHECK_EQUAL(6.7, m2(1, 0));
122}
123} // namespace StandardMatrixUnitTests
124
125namespace StandardMatrixOperationsUnitTests
126{
128{
129 double lhs_buf[] = {2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0,
130 18.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 32.0};
131 double rhs_buf[] = {3.0, 6.0, 9.0, 12.0, 15.0, 18.0, 21.0, 24.0,
132 27.0, 30.0, 33.0, 36.0, 39.0, 42.0, 45.0, 48.0};
133
134 NekMatrix<double> lhs1(4, 4, lhs_buf);
135 std::shared_ptr<NekMatrix<NekMatrix<double>, ScaledMatrixTag>> lhs2;
136 std::shared_ptr<NekMatrix<NekMatrix<double>, BlockMatrixTag>> lhs3;
137
138 NekMatrix<double> rhs1(4, 4, rhs_buf);
139 std::shared_ptr<NekMatrix<NekMatrix<double>, ScaledMatrixTag>> rhs2;
140 std::shared_ptr<NekMatrix<NekMatrix<double>, BlockMatrixTag>> rhs3;
141
142 GenerateMatrices(lhs1, 2.0, 2, 2, lhs2, lhs3);
143 GenerateMatrices(rhs1, 3.0, 2, 2, rhs2, rhs3);
144
145 double result_buf[] = {5, 10, 15, 20, 25, 30, 35, 40,
146 45, 50, 55, 60, 65, 70, 75, 80};
147 NekMatrix<double> result(4, 4, result_buf);
148
149 RunAllAddCombinations(lhs1, *lhs2, *lhs3, rhs1, *rhs2, *rhs3, result);
150
151 NekMatrix<double> rhs_transposed = Transpose(rhs1);
152 NekMatrix<double> result_with_one_operand_transposed =
153 lhs1 + rhs_transposed;
154
155 double expected_result_with_one_operand_transposed[] = {
156 5, 19, 33, 47, 16, 30, 44, 58, 27, 41, 55, 69, 38, 52, 66, 80};
157 BOOST_CHECK_EQUAL(
158 NekMatrix<double>(4, 4, expected_result_with_one_operand_transposed),
159 result_with_one_operand_transposed);
160}
161
162BOOST_AUTO_TEST_CASE(TestSubtraction)
163{
164 double lhs_buf[] = {2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0,
165 18.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 32.0};
166 double rhs_buf[] = {3.0, 6.0, 9.0, 12.0, 15.0, 18.0, 21.0, 24.0,
167 27.0, 30.0, 33.0, 36.0, 39.0, 42.0, 45.0, 48.0};
168
169 NekMatrix<double> lhs1(4, 4, lhs_buf);
170 std::shared_ptr<NekMatrix<NekMatrix<double>, ScaledMatrixTag>> lhs2;
171 std::shared_ptr<NekMatrix<NekMatrix<double>, BlockMatrixTag>> lhs3;
172
173 NekMatrix<double> rhs1(4, 4, rhs_buf);
174 std::shared_ptr<NekMatrix<NekMatrix<double>, ScaledMatrixTag>> rhs2;
175 std::shared_ptr<NekMatrix<NekMatrix<double>, BlockMatrixTag>> rhs3;
176
177 GenerateMatrices(lhs1, 2.0, 2, 2, lhs2, lhs3);
178 GenerateMatrices(rhs1, 3.0, 2, 2, rhs2, rhs3);
179
180 double result_buf[] = {-1, -2, -3, -4, -5, -6, -7, -8,
181 -9, -10, -11, -12, -13, -14, -15, -16};
182 NekMatrix<double> result(4, 4, result_buf);
183
184 RunAllSubCombinations(lhs1, *lhs2, *lhs3, rhs1, *rhs2, *rhs3, result);
185}
186
187BOOST_AUTO_TEST_CASE(TestThreeAdditions)
188{
189 double lhs_buf[] = {2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0,
190 18.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 32.0};
191 double middle_buf[] = {2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0,
192 18.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 32.0};
193 double rhs_buf[] = {3.0, 6.0, 9.0, 12.0, 15.0, 18.0, 21.0, 24.0,
194 27.0, 30.0, 33.0, 36.0, 39.0, 42.0, 45.0, 48.0};
195
196 NekMatrix<double> lhs(4, 4, lhs_buf);
197 NekMatrix<double> middle(4, 4, middle_buf);
198 NekMatrix<double> rhs(4, 4, rhs_buf);
199
200 NekMatrix<double> result = lhs + middle + rhs;
201}
202
203BOOST_AUTO_TEST_CASE(TestTransposedMultiplication)
204{
205 {
206 double buf1[] = {1, 4, 2, 5, 3, 6};
207 double buf2[] = {10, 11, 12, 13, 14, 15};
208
209 double transposed_buf1[] = {1, 2, 3, 4, 5, 6};
210
211 double transposed_buf2[] = {10, 13, 11, 14, 12, 15};
212
213 double result_buf[] = {68, 167, 86, 212};
214
215 double transposed_result_buf[] = {62, 85, 108, 67, 92,
216 117, 72, 99, 126};
217
218 NekMatrix<double> m1(2, 3, buf1);
219 NekMatrix<double> m2(3, 2, buf2);
220 NekMatrix<double> transposed_m1(3, 2, transposed_buf1);
221 NekMatrix<double> transposed_m2(2, 3, transposed_buf2);
222
223 NekMatrix<double> expected_result(2, 2, result_buf);
224 NekMatrix<double> expected_transposed_result(3, 3,
225 transposed_result_buf);
226
227 BOOST_CHECK_EQUAL(expected_result, m1 * m2);
228 BOOST_CHECK_EQUAL(expected_transposed_result,
229 transposed_m1 * transposed_m2);
230 m1.Transpose();
231 m2.Transpose();
232 BOOST_CHECK_EQUAL(expected_transposed_result, m1 * m2);
233 m1.Transpose();
234 m2.Transpose();
235 BOOST_CHECK_EQUAL(expected_result, m1 * m2);
236
239 BOOST_CHECK_EQUAL(expected_result, m1 * m2);
240 BOOST_CHECK_EQUAL(expected_transposed_result, tm1 * tm2);
241 }
242}
243
244BOOST_AUTO_TEST_CASE(TestGlobalTransposeMethod)
245{
246 double buf1[] = {1, 4, 2, 5, 3, 6};
247 double buf2[] = {10, 11, 12, 13, 14, 15};
248
249 double transposed_buf1[] = {1, 2, 3, 4, 5, 6};
250
251 double transposed_buf2[] = {10, 13, 11, 14, 12, 15};
252
253 NekMatrix<double> m1(2, 3, buf1);
254 NekMatrix<double> m2(3, 2, buf2);
255
258
259 NekMatrix<double> expected_tm1(3, 2, transposed_buf1);
260 NekMatrix<double> expected_tm2(2, 3, transposed_buf2);
261
262 BOOST_CHECK_EQUAL(tm1, expected_tm1);
263 BOOST_CHECK_EQUAL(tm2, expected_tm2);
264
265 BOOST_CHECK_EQUAL(m1(0, 1), 2);
266 BOOST_CHECK_EQUAL(tm1(1, 0), 2);
267
268 m1(0, 1) = 89;
269 BOOST_CHECK_EQUAL(m1(0, 1), 89);
270 BOOST_CHECK_EQUAL(tm1(1, 0), 89);
271
272 BOOST_CHECK_EQUAL(m1.GetTransposeFlag(), 'N');
273 BOOST_CHECK_EQUAL(tm1.GetTransposeFlag(), 'T');
274}
275
276BOOST_AUTO_TEST_CASE(TestScalarMatrixMultiply)
277{
278 double buf1[] = {1, 4, 2, 5, 3, 6};
279
280 NekMatrix<double> m1(2, 3, buf1);
281 NekMatrix<double> m2 = m1 * 2.0;
282 NekMatrix<double> m3 = 3.0 * m1;
283
284 double m2_expected_result_buf[] = {2.0, 8.0, 4.0, 10.0, 6.0, 12.0};
285
286 double m3_expected_result_buf[] = {3.0, 12.0, 6.0, 15.0, 9.0, 18.0};
287
288 NekMatrix<double> m2_expected_result(2, 3, m2_expected_result_buf);
289 NekMatrix<double> m3_expected_result(2, 3, m3_expected_result_buf);
290
291 BOOST_CHECK_EQUAL(m2_expected_result, m2);
292 BOOST_CHECK_EQUAL(m3_expected_result, m3);
293}
294
295} // namespace StandardMatrixOperationsUnitTests
296
297} // namespace Nektar
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)
void RunAllAddCombinations(const NekMatrix< NekDouble, StandardMatrixTag > &l1, const NekMatrix< NekMatrix< NekDouble, LhsScaledInnerMatrixType >, ScaledMatrixTag > &l2, const NekMatrix< NekMatrix< NekDouble, LhsBlockInnerMatrixType >, BlockMatrixTag > &l3, const NekMatrix< NekDouble, StandardMatrixTag > &r1, const NekMatrix< NekMatrix< NekDouble, RhsScaledInnerMatrixType >, ScaledMatrixTag > &r2, const NekMatrix< NekMatrix< NekDouble, RhsBlockInnerMatrixType >, BlockMatrixTag > &r3, const NekMatrix< NekDouble, StandardMatrixTag > &result)
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
void RunAllSubCombinations(const NekMatrix< NekDouble, StandardMatrixTag > &l1, const NekMatrix< NekMatrix< NekDouble, LhsScaledInnerMatrixType >, ScaledMatrixTag > &l2, const NekMatrix< NekMatrix< NekDouble, LhsBlockInnerMatrixType >, BlockMatrixTag > &l3, const NekMatrix< NekDouble, StandardMatrixTag > &r1, const NekMatrix< NekMatrix< NekDouble, RhsScaledInnerMatrixType >, ScaledMatrixTag > &r2, const NekMatrix< NekMatrix< NekDouble, RhsBlockInnerMatrixType >, BlockMatrixTag > &r3, const NekMatrix< NekDouble, StandardMatrixTag > &result)