Nektar++
testLinearSystem.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: testLinearSystem.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: Test code for NekVector
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <LibUtilities/LinearAlgebra/NekLinSys.hpp>
36#include <boost/test/tools/floating_point_comparison.hpp>
37#include <boost/test/unit_test.hpp>
38
39namespace Nektar
40{
41namespace LinearSystemUnitTests
42{
43BOOST_AUTO_TEST_CASE(TestLinearSystemSolveWithReturnValue)
44{
45 double matrix_buf[] = {10, 5, 2};
46
47 double b_buf[] = {20, 50, 10};
48
49 std::shared_ptr<NekMatrix<double>> A(
50 new NekMatrix<double>(3, 3, matrix_buf, eDIAGONAL));
51 NekVector<double> b1(3, b_buf);
52 std::shared_ptr<NekVector<double>> b2(new NekVector<double>(3, b_buf));
53 NekVector<double> b3(3, b_buf);
54 std::shared_ptr<NekVector<double>> b4(new NekVector<double>(3, b_buf));
55
56 LinearSystem linsys(A);
57
58 NekVector<double> result1 = linsys.Solve(b1);
59 NekVector<double> result2 = linsys.Solve(b2);
60 NekVector<double> result3 = linsys.Solve(&b1);
61 NekVector<double> result4 = linsys.Solve(b2.get());
62
63 NekVector<double> result5 = linsys.Solve(b3);
64 NekVector<double> result6 = linsys.Solve(b4);
65 NekVector<double> result7 = linsys.Solve(&b3);
66 NekVector<double> result8 = linsys.Solve(b4.get());
67
68 double expected_result_buf[] = {2, 10, 5};
69 NekVector<double> expectedResult(3, expected_result_buf);
70
71 BOOST_CHECK_EQUAL(result1, expectedResult);
72 BOOST_CHECK_EQUAL(result2, expectedResult);
73 BOOST_CHECK_EQUAL(result3, expectedResult);
74 BOOST_CHECK_EQUAL(result4, expectedResult);
75 BOOST_CHECK_EQUAL(result5, expectedResult);
76 BOOST_CHECK_EQUAL(result6, expectedResult);
77 BOOST_CHECK_EQUAL(result7, expectedResult);
78 BOOST_CHECK_EQUAL(result8, expectedResult);
79}
80
81BOOST_AUTO_TEST_CASE(TestLinearSystemSolveWithSoluationAsParameter)
82{
83 double matrix_buf[] = {10, 5, 2};
84
85 double b_buf[] = {20, 50, 10};
86
87 std::shared_ptr<NekMatrix<double>> A(
88 new NekMatrix<double>(3, 3, matrix_buf, eDIAGONAL));
89 NekVector<double> b1(3, b_buf);
90 std::shared_ptr<NekVector<double>> b2(new NekVector<double>(3, b_buf));
91 NekVector<double> b3(3, b_buf);
92 std::shared_ptr<NekVector<double>> b4(new NekVector<double>(3, b_buf));
93
94 NekVector<double> result1(3, 0.0);
95 std::shared_ptr<NekVector<double>> result2(new NekVector<double>(3, 0.0));
96 NekVector<double> result3(3, 0.0);
97 std::shared_ptr<NekVector<double>> result4(new NekVector<double>(3, 0.0));
98
99 LinearSystem linsys(A);
100 double expected_result_buf[] = {2, 10, 5};
101 NekVector<double> expectedResult(3, expected_result_buf);
102
103 linsys.Solve(b1, result1);
104 linsys.Solve(b2, result3);
105 linsys.Solve(b3, result4);
106 linsys.Solve(b4, result2);
107
108 BOOST_CHECK_EQUAL(result1, expectedResult);
109 BOOST_CHECK_EQUAL(*result2, expectedResult);
110 BOOST_CHECK_EQUAL(result3, expectedResult);
111 BOOST_CHECK_EQUAL(*result4, expectedResult);
112}
113
114BOOST_AUTO_TEST_CASE(TestDiagonalSystem)
115{
116 double matrix_buf[] = {10, 5, 2};
117
118 double result_buf[] = {20, 50, 10};
119
120 std::shared_ptr<NekMatrix<double>> A(
121 new NekMatrix<double>(3, 3, matrix_buf, eDIAGONAL));
122 std::shared_ptr<NekVector<double>> b(new NekVector<double>(3, result_buf));
123
124 LinearSystem linsys(A);
125
126 NekVector<double> result = linsys.Solve(b);
127
128 double expected_result_buf[] = {2, 10, 5};
129 NekVector<double> expectedResult(3, expected_result_buf);
130
131 BOOST_CHECK_EQUAL(result, expectedResult);
132 std::shared_ptr<NekVector<double>> b1 = b;
133 NekVector<double> result1 = linsys.Solve(b1);
134 BOOST_CHECK_EQUAL(result1, expectedResult);
135
136 NekVector<double> result2 = linsys.SolveTranspose(b);
137 BOOST_CHECK_EQUAL(expectedResult, result2);
138}
139
140BOOST_AUTO_TEST_CASE(TestDiagonalSystemAsFullSystem)
141{
142 double matrix_buf[] = {10.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 2.0};
143
144 double result_buf[] = {20, 50, 10};
145
146 std::shared_ptr<NekMatrix<double>> A(
147 new NekMatrix<double>(3, 3, matrix_buf, eFULL));
148 std::shared_ptr<NekVector<double>> b(new NekVector<double>(3, result_buf));
149
150 LinearSystem linsys(A);
151
152 NekVector<double> result = linsys.Solve(b);
153
154 double epsilon = 1e-14;
155
156 BOOST_CHECK_CLOSE(result[0], 2.0, epsilon);
157 BOOST_CHECK_CLOSE(result[1], 10.0, epsilon);
158 BOOST_CHECK_CLOSE(result[2], 5.0, epsilon);
159}
160
161BOOST_AUTO_TEST_CASE(TestSmallFullSystem)
162{
163
164 double matrix_buf[] = {81, -5, -28, 4};
165
166 double b_buf[] = {-941, 348};
167
168 std::shared_ptr<NekMatrix<double>> A(
169 new NekMatrix<double>(2, 2, matrix_buf, eFULL));
170 std::shared_ptr<NekVector<double>> b(new NekVector<double>(2, b_buf));
171 LinearSystem linsys(A);
172
173 NekVector<double> result = linsys.Solve(b);
174 double epsilon = 1e-11;
175
176 BOOST_CHECK_CLOSE(result[0], 32.5, epsilon);
177 BOOST_CHECK_CLOSE(result[1], 127.625, epsilon);
178
179 result = linsys.SolveTranspose(b);
180 BOOST_CHECK_CLOSE(result[0], -11.0, epsilon);
181 BOOST_CHECK_CLOSE(result[1], 10.0, epsilon);
182}
183
184BOOST_AUTO_TEST_CASE(TestLargeFullSystem)
185{
186
187 // Larger matrix.
188 double matrix_buf[] = {
189 -85, -55, -37, -35, 97, 50, 79, 56, 49, 63, 57, -59, 45,
190 -8, -93, 92, 43, -62, 77, 66, 54, -5, 99, -61, -50, -12,
191 -18, 31, -26, -62, 1, -47, -91, -47, -61, 41, -58, -90, 53,
192 -1, 94, 83, -86, 23, -84, 19, -50, 88, -53, 85, 49, 78,
193 17, 72, -99, -85, -86, 30, 80, 72, 66, -29, -91, -53, -19,
194 -47, 68, -72, -87, 79, 43, -66, -53, -61, -23, -37, 31, -34,
195 -42, 88, -76, -65, 25, 28, -61, -60, 9, 29, -66, -32, 78,
196 39, 94, 68, -17, -98, -36, 40, 22, 5};
197
198 double b_buf[] = {12719, -3169, -16810, 7408, -14945,
199 -6822, 10166, 7023, 8679, -11826};
200
201 std::shared_ptr<NekMatrix<double>> A(
202 new NekMatrix<double>(10, 10, matrix_buf, eFULL));
203 std::shared_ptr<NekVector<double>> b(new NekVector<double>(10, b_buf));
204 LinearSystem linsys(A);
205
206 NekVector<double> result = linsys.SolveTranspose(b);
207 double epsilon = 1e-11;
208 BOOST_CHECK_CLOSE(result[0], -88.0, epsilon);
209 BOOST_CHECK_CLOSE(result[1], -43.0, epsilon);
210 BOOST_CHECK_CLOSE(result[2], -73.0, epsilon);
211 BOOST_CHECK_CLOSE(result[3], 25.0, epsilon);
212 BOOST_CHECK_CLOSE(result[4], 4.0, epsilon);
213 BOOST_CHECK_CLOSE(result[5], -59.0, epsilon);
214 BOOST_CHECK_CLOSE(result[6], 62.0, epsilon);
215 BOOST_CHECK_CLOSE(result[7], -55.0, epsilon);
216 BOOST_CHECK_CLOSE(result[8], 25.0, epsilon);
217 BOOST_CHECK_CLOSE(result[9], 9.0, epsilon);
218}
219
220// void testSolvingBlockDiagonalMatrices()
221// {
222// //typedef NekMatrix<double, DiagonalMatrixTag, eBlock>
223// BlockMatrix;
224// //typedef BlockMatrix::InnerMatrixType InnerMatrixType;
225// //
226// //std::shared_ptr<BlockMatrix> m1(new BlockMatrix(4, 2, 2));
227// //const double m1_buf0[] = {-91, -47, 94, 83};
228// //const double m1_buf1[] = {49, 78, 80, 72};
229// //const double m1_buf2[] = {87, 79, 31, -34};
230// //const double m1_buf3[] = {9, 29, -17, -98};
231// //
232// //m1->GetBlock(0,0) = InnerMatrixType(2, 2, m1_buf0);
233// //m1->GetBlock(1,1) = InnerMatrixType(2, 2, m1_buf1);
234// //m1->GetBlock(2,2) = InnerMatrixType(2, 2, m1_buf2);
235// //m1->GetBlock(3,3) = InnerMatrixType(2, 2, m1_buf3);
236// //
237// //double b_vals[] = {-7198, 7642, 1344, 3744, -9968, 115, -2469,
238// 8424};
239// //std::shared_ptr<NekVector<double> > b(new NekVector<double>(8,
240// b_vals));
241// //
242// //LinearSystem<BlockMatrix, NekVector<double> > linsys(m1, b);
243// //
244// //NekVector<double> result = linsys.Solve();
245// //double epsilon = 1e-11;
246// //BOOST_CHECK_CLOSE(result[0], 76.0, epsilon);
247// //BOOST_CHECK_CLOSE(result[1], 6.0, epsilon);
248// //BOOST_CHECK_CLOSE(result[2], 72.0, epsilon);
249// //BOOST_CHECK_CLOSE(result[3], -28.0, epsilon);
250// //BOOST_CHECK_CLOSE(result[4], -61.0, epsilon);
251// //BOOST_CHECK_CLOSE(result[5], -59.0, epsilon);
252// //BOOST_CHECK_CLOSE(result[6], 6.0, epsilon);
253// //BOOST_CHECK_CLOSE(result[7], -87.0, epsilon);
254// }
255
256BOOST_AUTO_TEST_CASE(TestFullSystemWithWrappedVectors)
257{
258 // Larger matrix.
259 double matrix_buf[] = {
260 -85, -55, -37, -35, 97, 50, 79, 56, 49, 63, 57, -59, 45,
261 -8, -93, 92, 43, -62, 77, 66, 54, -5, 99, -61, -50, -12,
262 -18, 31, -26, -62, 1, -47, -91, -47, -61, 41, -58, -90, 53,
263 -1, 94, 83, -86, 23, -84, 19, -50, 88, -53, 85, 49, 78,
264 17, 72, -99, -85, -86, 30, 80, 72, 66, -29, -91, -53, -19,
265 -47, 68, -72, -87, 79, 43, -66, -53, -61, -23, -37, 31, -34,
266 -42, 88, -76, -65, 25, 28, -61, -60, 9, 29, -66, -32, 78,
267 39, 94, 68, -17, -98, -36, 40, 22, 5};
268
269 double b_buf[] = {12719, -3169, -16810, 7408, -14945,
270 -6822, 10166, 7023, 8679, -11826};
271 double result_buf[10];
272
273 std::shared_ptr<NekMatrix<double>> A(
274 new NekMatrix<double>(10, 10, matrix_buf, eFULL));
275 std::shared_ptr<NekVector<double>> b(
276 new NekVector<double>(Array<OneD, double>(10, b_buf), eWrapper));
277 NekVector<double> result(Array<OneD, double>(10, result_buf), eWrapper);
278
279 LinearSystem linsys(A);
280
281 linsys.SolveTranspose(b, result);
282 double epsilon = 1e-11;
283 BOOST_CHECK_CLOSE(result[0], -88.0, epsilon);
284 BOOST_CHECK_CLOSE(result[1], -43.0, epsilon);
285 BOOST_CHECK_CLOSE(result[2], -73.0, epsilon);
286 BOOST_CHECK_CLOSE(result[3], 25.0, epsilon);
287 BOOST_CHECK_CLOSE(result[4], 4.0, epsilon);
288 BOOST_CHECK_CLOSE(result[5], -59.0, epsilon);
289 BOOST_CHECK_CLOSE(result[6], 62.0, epsilon);
290 BOOST_CHECK_CLOSE(result[7], -55.0, epsilon);
291 BOOST_CHECK_CLOSE(result[8], 25.0, epsilon);
292 BOOST_CHECK_CLOSE(result[9], 9.0, epsilon);
293}
294
295BOOST_AUTO_TEST_CASE(TestSymmetricMatrix)
296{
297 NekMatrix<double, StandardMatrixTag> m(3, 3, eSYMMETRIC);
298 m.SetValue(0, 0, 3.0);
299 m.SetValue(0, 1, 4.0);
300 m.SetValue(0, 2, 9);
301 m.SetValue(1, 1, 7);
302 m.SetValue(1, 2, 2);
303 m.SetValue(2, 2, 11);
304
305 BOOST_CHECK_EQUAL(4.0, m(1, 0));
306
307 double b_buf[] = {-20.0, 30.0, -20.0};
308 NekVector<double> b(3, b_buf);
309
310 LinearSystem linsys(m);
311 NekVector<double> x = linsys.Solve(b);
312 NekVector<double> xt = linsys.SolveTranspose(b);
313
314 double epsilon = 1e-11;
315 BOOST_CHECK_CLOSE(x[0], 3.0, epsilon);
316 BOOST_CHECK_CLOSE(x[1], 4.0, epsilon);
317 BOOST_CHECK_CLOSE(x[2], -5.0, epsilon);
318 BOOST_CHECK_CLOSE(xt[0], 3.0, epsilon);
319 BOOST_CHECK_CLOSE(xt[1], 4.0, epsilon);
320 BOOST_CHECK_CLOSE(xt[2], -5.0, epsilon);
321}
322} // namespace LinearSystemUnitTests
323} // namespace Nektar
Nektar::LinearSystem::SolveTranspose
RawType_t< VectorType > SolveTranspose(const VectorType &b)
Definition: NekLinSys.hpp:475
Nektar::LinearSystem::Solve
RawType_t< VectorType > Solve(const VectorType &b)
Definition: NekLinSys.hpp:452
Nektar::LinearSystemUnitTests::BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(TestLinearSystemSolveWithReturnValue)
Definition: testLinearSystem.cpp:43
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
Generated on Wed Jan 31 2024 08:26:49 for Nektar++ by doxygen 1.9.4