1 /* Siconos is a program dedicated to modeling, simulation and control
2 * of non smooth dynamical systems.
3 *
4 * Copyright 2021 INRIA.
5 *
6 * Licensed under the Apache License, Version 2.0 (the "License");
7 * you may not use this file except in compliance with the License.
8 * You may obtain a copy of the License at
9 *
10 * http://www.apache.org/licenses/LICENSE-2.0
11 *
12 * Unless required by applicable law or agreed to in writing, software
13 * distributed under the License is distributed on an "AS IS" BASIS,
14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15 * See the License for the specific language governing permissions and
16 * limitations under the License.
17 */
18 #include "SiconosConfig.h"
19
20 #include "SimpleMatrixTest.hpp"
21 #include "SiconosAlgebra.hpp"
22 #include "SiconosAlgebraProd.hpp"
23 #include "SiconosAlgebraScal.hpp"
24 #include "SimpleMatrixFriends.hpp"
25 #include "SiconosMatrixSetBlock.hpp"
26 #include "SiconosVectorFriends.hpp"
27 #include "NumericsMatrix.h"
28 #include "NumericsSparseMatrix.h"
29
30 #include <boost/numeric/ublas/matrix_sparse.hpp>
31
32 #define CPPUNIT_ASSERT_NOT_EQUAL(message, alpha, omega) \
33 if ((alpha) == (omega)) CPPUNIT_FAIL(message);
34
35 CPPUNIT_TEST_SUITE_REGISTRATION(SimpleMatrixTest);
36
37 using namespace Siconos;
38
39 #define DEBUG_MESSAGES
40 #include "siconos_debug.h"
41
42
43 // Note FP: tests are (rather) complete for Dense objects but many are missing for other cases (Triang, Symm etc ...).
44
setUp()45 void SimpleMatrixTest::setUp()
46 {
47 tol = 1e-9;
48
49 fic1 = "mat1.dat"; // 2 X 2
50 fic2 = "mat2.dat"; // 2 X 3
51 SicM.reset(new SimpleMatrix(fic1, 1));
52 SimM.reset(new SimpleMatrix(fic2, 1));
53
54 std::vector<double> v3(2, 0);
55 std::vector<double> v4(2, 0);
56 std::vector<double> v5(3, 0);
57 v4[0] = 6;
58 v4[1] = 9;
59 v5[0] = 8;
60 v5[1] = 9;
61 v5[2] = 10;
62
63
64
65 vect1.reset(new SiconosVector(v3));
66 vect2.reset(new SiconosVector(v4)); // vect2 != vect1, but vect2 == SimM second column
67 vect3.reset(new SiconosVector(v5)); // vect3 != vect1, but vect3 == SimM second row
68
69 // Dense
70 D.reset(new DenseMat(2, 2));
71 for(unsigned i = 0; i < D->size1(); ++ i)
72 for(unsigned j = 0; j < D->size2(); ++ j)
73 (*D)(i, j) = 3 * i + j;
74
75 // Triang
76 T.reset(new TriangMat(3, 3));
77 for(unsigned i = 0; i < T->size1(); ++ i)
78 for(unsigned j = i; j < T->size2(); ++ j)
79 (*T)(i, j) = 3 * i + j;
80 T2.reset(new TriangMat(4, 4));
81 for(unsigned i = 0; i < T2->size1(); ++ i)
82 for(unsigned j = i; j < T2->size2(); ++ j)
83 (*T2)(i, j) = 3 * i + j;
84
85 // Sym
86 S.reset(new SymMat(3, 3));
87 for(unsigned i = 0; i < S->size1(); ++ i)
88 for(unsigned j = i; j < S->size2(); ++ j)
89 (*S)(i, j) = 3 * i + j;
90 S2.reset(new SymMat(4, 4));
91 for(unsigned i = 0; i < S2->size1(); ++ i)
92 for(unsigned j = i; j < S2->size2(); ++ j)
93 (*S2)(i, j) = 3 * i + j;
94
95 // Sparse
96 SP.reset(new SparseMat(4, 4));
97 for(unsigned i = 0; i < SP->size1(); ++ i)
98 for(unsigned j = 0; j < SP->size2(); ++ j)
99 (*SP)(i, j) = 3 * i + j;
100
101 SP2.reset(new SparseMat(4, 4));
102 for(unsigned i = 0; i < SP2->size1(); ++ i)
103 for(unsigned j = 0; j < SP->size2()-1; ++ j)
104 if(i != j)
105 (*SP2)(i, j) = 3 * i + j;
106
107 SP3.reset(new SparseMat(3, 3));
108 (*SP3)(0,0) = 1.0;
109 (*SP3)(0,2) = 2.0;
110 (*SP3)(1,2) = 3.0;
111 (*SP3)(2,0) = 4.0;
112 (*SP3)(2,1) = 5.0;
113 (*SP3)(2,2) = 5.0;
114
115 SP4.reset(new SparseMat(3, 3));
116 (*SP4)(0,0) = 1.0;
117 (*SP4)(0,2) = 4.0;
118 (*SP4)(1,1) = 1.0;
119 (*SP4)(1,2) = 5.0;
120 (*SP4)(2,0) = 4.0;
121 (*SP4)(2,1) = 5.0;
122 (*SP4)(2,2) = 77.0;
123
124
125 // Sparse Coordinate
126 SP_coor.reset(new SparseCoordinateMat(4, 4));
127 for(unsigned i = 0; i < SP->size1(); ++ i)
128 for(unsigned j = 0; j < SP->size2(); ++ j)
129 (*SP_coor)(i, j) = 3 * i + j;
130
131 // Banded
132 Band.reset(new BandedMat(4, 4, 1, 1));
133 for(signed i = 0; i < signed(Band->size1()); ++ i)
134 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
135 (*Band)(i, j) = 3 * i + j;
136 Band2.reset(new BandedMat(4, 3, 1, 1));
137 for(signed i = 0; i < signed(Band2->size1()); ++ i)
138 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band2->size2())); ++ j)
139 (*Band2)(i, j) = 3 * i + j;
140
141 // Zero
142 Z.reset(new ZeroMat(3, 3));
143 Z2.reset(new ZeroMat(4, 4));
144 // Identity
145 I.reset(new IdentityMat(3, 3));
146 I2.reset(new IdentityMat(4, 4));
147
148 // BlockMat
149 size = 10;
150 size2 = 10;
151
152 C.reset(new SimpleMatrix(size, size));
153 A.reset(new SimpleMatrix("A.dat"));
154 B.reset(new SimpleMatrix("B.dat"));
155
156 m1.reset(new SimpleMatrix(size - 2, size - 2, 1));
157 m2.reset(new SimpleMatrix(size - 2, 2, 2));
158 m3.reset(new SimpleMatrix(2, size - 2, 3));
159 m4.reset(new SimpleMatrix(2, 2, 4));
160 m5.reset(new SimpleMatrix(size2 - 2, size2 - 2, 1));
161 m6.reset(new SimpleMatrix(size2 - 2, 2, 2));
162 m7.reset(new SimpleMatrix(2, size2 - 2, 3));
163 m8.reset(new SimpleMatrix(2, 2, 4));
164 Ab.reset(new BlockMatrix(m1, m2, m3, m4));
165 Bb.reset(new BlockMatrix(3 * *Ab));
166 Cb.reset(new BlockMatrix(m5, m6, m7, m8));
167
168
169 }
170
tearDown()171 void SimpleMatrixTest::tearDown()
172 {}
173
174 //______________________________________________________________________________
175
testConstructor0()176 void SimpleMatrixTest::testConstructor0() // constructor with TYP and dim
177 {
178 std::cout << "====================================" <<std::endl;
179 std::cout << "=== Simple Matrix tests start ...=== " <<std::endl;
180 std::cout << "====================================" <<std::endl;
181 std::cout << "--> Test: constructor 0." <<std::endl;
182 SP::SimpleMatrix test(new SimpleMatrix(2, 3));
183
184 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor0 : ", test->num() == Siconos::DENSE, true);
185 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor0 : ", test->size(0) == 2, true);
186 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor0 : ", test->size(1) == 3, true);
187 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor0 : ", test->normInf() < tol, true);
188 std::cout << "--> Constructor 0 test ended with success." <<std::endl;
189 }
190
testConstructor1()191 void SimpleMatrixTest::testConstructor1() // Copy constructor, from a SimpleMatrix
192 {
193 std::cout << "--> Test: constructor 1." <<std::endl;
194 SP::SimpleMatrix test(new SimpleMatrix(*SimM));
195 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor1 : ", *test == *SimM, true);
196 std::cout << "--> Constructor 1 (copy) test ended with success." <<std::endl;
197 }
198
testConstructor2()199 void SimpleMatrixTest::testConstructor2() // Copy constructor, from a SiconosMatrix
200 {
201 std::cout << "--> Test: constructor 2." <<std::endl;
202 SP::SimpleMatrix test(new SimpleMatrix(*SicM));
203 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor2 : ", *test == *SicM, true);
204 std::cout << "--> Constructor 2 (copy) test ended with success." <<std::endl;
205 }
206
testConstructor3()207 void SimpleMatrixTest::testConstructor3() // Copy constructor, from a BlockMatrix
208 {
209 std::cout << "--> Test: constructor 3." <<std::endl;
210 SP::SimpleMatrix test(new SimpleMatrix(*Ab));
211 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor3 : ", *test == *Ab, true);
212 std::cout << "--> Constructor 3 (copy) test ended with success." <<std::endl;
213 }
214
215
testConstructor4()216 void SimpleMatrixTest::testConstructor4()
217 {
218 std::cout << "--> Test: constructor 4." <<std::endl;
219 SP::SiconosMatrix test(new SimpleMatrix(*D));
220 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor4 : ", test->num() == Siconos::DENSE, true);
221 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor4 : ", norm_inf(test->getDense() - *D) == 0, true);
222 std::cout << "--> Constructor 4 test ended with success." <<std::endl;
223 }
224
testConstructor5()225 void SimpleMatrixTest::testConstructor5()
226 {
227 std::cout << "--> Test: constructor 5." <<std::endl;
228 SP::SiconosMatrix test(new SimpleMatrix(*T));
229 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor5 : ", test->num() == Siconos::TRIANGULAR, true);
230 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor5 : ", norm_inf(test->getTriang() - *T) == 0, true);
231 std::cout << "--> Constructor 5 test ended with success." <<std::endl;
232 }
233
testConstructor6()234 void SimpleMatrixTest::testConstructor6()
235 {
236 std::cout << "--> Test: constructor 6." <<std::endl;
237 SP::SiconosMatrix test(new SimpleMatrix(*S));
238 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor6 : ", test->num() == Siconos::SYMMETRIC, true);
239 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor6 : ", norm_inf(test->getSym() - *S) == 0, true);
240 std::cout << "--> Constructor 6 test ended with success." <<std::endl;
241 }
242
testConstructor7()243 void SimpleMatrixTest::testConstructor7()
244 {
245 std::cout << "--> Test: constructor 7." <<std::endl;
246 SP::SiconosMatrix test(new SimpleMatrix(*SP));
247 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor7 : ", test->num() == Siconos::SPARSE, true);
248 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor7 : ", norm_inf(test->getSparse() - *SP) == 0, true);
249 std::cout << "--> Constructor 7 test ended with success." <<std::endl;
250 }
251
testConstructor8()252 void SimpleMatrixTest::testConstructor8()
253 {
254 std::cout << "--> Test: constructor 8." <<std::endl;
255 std::cout << "--> Constructor 8 test ended with success." <<std::endl;
256 SP::SiconosMatrix test(new SimpleMatrix(*Band));
257 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor8 : ", test->num() == Siconos::BANDED, true);
258 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor8 : ", norm_inf(test->getBanded() - *Band) == 0, true);
259 }
260
testConstructor9()261 void SimpleMatrixTest::testConstructor9() // constructor with TYP and dim and input value
262 {
263 std::cout << "--> Test: constructor 9." <<std::endl;
264 SP::SimpleMatrix test(new SimpleMatrix(2, 3, 4.5));
265
266 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor9 : ", test->num() == Siconos::DENSE, true);
267 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor9 : ", test->size(0) == 2, true);
268 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor9 : ", test->size(1) == 3, true);
269 for(unsigned int i = 0; i < 2; ++i)
270 for(unsigned int j = 0 ; j < 3; ++j)
271 {
272 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor9 : ", (*test)(i, j) == 4.5, true);
273 }
274 std::cout << "--> Constructor 9 test ended with success." <<std::endl;
275 }
276
testConstructor10()277 void SimpleMatrixTest::testConstructor10()
278 {
279 std::cout << "--> Test: constructor 10." <<std::endl;
280 SP::SiconosMatrix test(new SimpleMatrix(fic1));
281 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor10 : ", *test == *SicM, true);
282 std::cout << "--> Constructor 10 test ended with success." <<std::endl;
283 }
284
testConstructor11()285 void SimpleMatrixTest::testConstructor11()
286 {
287 std::cout << "--> Test: constructor 11." <<std::endl;
288 std::cout << "--> Constructor 11 test ended with success." <<std::endl;
289 SP::SiconosMatrix test(new SimpleMatrix(*Z));
290 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor11 : ", test->num() == Siconos::ZERO, true);
291 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor11 : ", test->normInf() == 0, true);
292 }
293
testConstructor12()294 void SimpleMatrixTest::testConstructor12()
295 {
296 std::cout << "--> Test: constructor 12." <<std::endl;
297 std::cout << "--> Constructor 12 test ended with success." <<std::endl;
298 SP::SiconosMatrix test(new SimpleMatrix(*I));
299 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor12 : ", test->num() == Siconos::IDENTITY, true);
300 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor12 : ", test->normInf() == 1, true);
301 }
302
testConstructor13()303 void SimpleMatrixTest::testConstructor13()
304 {
305 std::cout << "--> Test: constructor 13." <<std::endl;
306 SP::SimpleMatrix test(new SimpleMatrix(4,4,Siconos::SPARSE));
307 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor3 : ",test->num() == Siconos::SPARSE, true);
308 std::cout << "--> Constructor 13 test ended with success." <<std::endl;
309 }
testConstructor14()310 void SimpleMatrixTest::testConstructor14()
311 {
312 std::cout << "--> Test: constructor 14." <<std::endl;
313 SP::SimpleMatrix test(new SimpleMatrix(4,4,Siconos::SPARSE_COORDINATE));
314 CPPUNIT_ASSERT_EQUAL_MESSAGE("testConstructor14 : ",test->num() == Siconos::SPARSE_COORDINATE, true);
315 std::cout << "--> Constructor 14 test ended with success." <<std::endl;
316 }
317 // Add tests with getDense ...
318
319 // Add tests with getDense ...
320
testZero()321 void SimpleMatrixTest::testZero()
322 {
323 std::cout << "--> Test: zero." <<std::endl;
324 SP::SiconosMatrix tmp(new SimpleMatrix(*SimM));
325 tmp->zero();
326 unsigned int n1 = tmp->size(0);
327 unsigned int n2 = tmp->size(1);
328 for(unsigned int i = 0; i < n1; ++i)
329 for(unsigned int j = 0; j < n2; ++j)
330 CPPUNIT_ASSERT_EQUAL_MESSAGE("testZero : ", (*tmp)(i, j) == 0, true);
331 std::cout << "--> zero test ended with success." <<std::endl;
332 }
333
testEye()334 void SimpleMatrixTest::testEye()
335 {
336 std::cout << "--> Test: eye." <<std::endl;
337 SP::SiconosMatrix tmp(new SimpleMatrix(*SimM));
338 tmp->eye();
339 unsigned int n1 = tmp->size(0);
340 unsigned int n2 = tmp->size(1);
341 for(unsigned int i = 0; i < n1; ++i)
342 for(unsigned int j = 0; j < n2; ++j)
343 if(i != j)
344 CPPUNIT_ASSERT_EQUAL_MESSAGE("testEye : ", (*tmp)(i, j) == 0, true);
345 else
346 CPPUNIT_ASSERT_EQUAL_MESSAGE("testEye : ", (*tmp)(i, j) == 1, true);
347 std::cout << "--> eye test ended with success." <<std::endl;
348 }
349
testResize()350 void SimpleMatrixTest::testResize()
351 {
352 std::cout << "--> Test: resize." <<std::endl;
353 SP::SiconosMatrix tmp(new SimpleMatrix(*SicM));
354 tmp->resize(3, 4);
355 unsigned int n1 = SicM->size(0);
356 unsigned int n2 = SicM->size(1);
357 CPPUNIT_ASSERT_EQUAL_MESSAGE("testResize : ", tmp->size(0) == 3, true);
358 CPPUNIT_ASSERT_EQUAL_MESSAGE("testResize : ", tmp->size(1) == 4, true);
359
360 for(unsigned int i = 0; i < n1; ++i)
361 for(unsigned int j = 0; j < n2; ++j)
362 CPPUNIT_ASSERT_EQUAL_MESSAGE("testResize : ", fabs((*tmp)(i, j) - (*SicM)(i, j)) < tol, true);
363 // for(unsigned int i = n1; i<3; ++i)
364 // for(unsigned int j=0;j<4;++j)
365 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testResize : ", fabs((*tmp)(i,j)) < tol, true);
366 // for(unsigned int j = n2; j<4; ++j)
367 // for(unsigned int i=0;i<3;++i)
368 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testResize : ", fabs((*tmp)(i,j)) < tol, true)
369 ;
370 // Check the effect of bool = false (ie preserve == false in boost resize)
371 // tmp->resize(6,8, false);
372 // tmp->display();
373 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testResize : ", tmp->size(0) == 6, true);
374 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testResize : ", tmp->size(1) == 8, true);
375 // for(unsigned int i = 0; i<6; ++i)
376 // for(unsigned int j=0;j<8;++j)
377 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testResize : ", (*tmp)(i,j) == 0 , true);
378 // // Reduction ...
379 // tmp->resize(1,2, false);
380 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testResize : ", tmp->size(0) == 1, true);
381 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testResize : ", tmp->size(1) == 2, true);
382 // for(unsigned int i = 0; i<1; ++i)
383 // for(unsigned int j=0;j<2;++j)
384 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testResize : ", (*tmp)(i,j) == 0 , true);
385 std::cout << "--> resize test ended with success." <<std::endl;
386 }
387
testNormInf()388 void SimpleMatrixTest::testNormInf()
389 {
390 std::cout << "--> Test: normInf." <<std::endl;
391 double n = SicM->normInf();
392 CPPUNIT_ASSERT_EQUAL_MESSAGE("testNormInf: ", n == 7, true);
393 std::cout << "--> normInf test ended with success." <<std::endl;
394 }
395
testSetBlock()396 void SimpleMatrixTest::testSetBlock()
397 {
398 std::cout << "--> Test: testSetBlock." <<std::endl;
399
400 // Copy of a sub-block of a Simple into a Simple
401 SP::SiconosMatrix MIn(new SimpleMatrix(10, 10));
402 for(unsigned int i = 0; i < 10; ++i)
403 for(unsigned int j = 0 ; j < 10; ++j)
404 (*MIn)(i, j) = i + j;
405
406 SP::SiconosMatrix MOut(new SimpleMatrix(5, 5));
407
408 Index subDim(2);
409 Index subPos(4);
410 subDim[0] = 2;
411 subDim[1] = 3;
412 subPos[0] = 1;
413 subPos[1] = 2;
414 subPos[2] = 1;
415 subPos[3] = 2;
416
417 setBlock(MIn, MOut, subDim, subPos);
418
419 for(unsigned int i = subPos[2]; i < subPos[2] + subDim[0]; ++i)
420 for(unsigned int j = subPos[3] ; j < subPos[3] + subDim[1]; ++j)
421 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSetBlock: ", fabs((*MOut)(i, j) - (*MIn)(i, j)) < tol, true);
422
423 // Copy of a sub-block of a Simple into a Block
424 Cb->zero();
425 setBlock(MIn, Cb, subDim, subPos);
426
427 for(unsigned int i = subPos[2]; i < subPos[2] + subDim[0]; ++i)
428 for(unsigned int j = subPos[3] ; j < subPos[3] + subDim[1]; ++j)
429 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSetBlock: ", fabs((*Cb)(i, j) - (*MIn)(i, j)) < tol, true);
430
431 // Copy of a sub-block of a Block into a Simple
432
433 MOut.reset(new SimpleMatrix(5, 5));
434 setBlock(Ab, MOut, subDim, subPos);
435
436 for(unsigned int i = subPos[2]; i < subPos[2] + subDim[0]; ++i)
437 for(unsigned int j = subPos[3] ; j < subPos[3] + subDim[1]; ++j)
438 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSetBlock: ", fabs((*MOut)(i, j) - (*Ab)(i, j)) < tol, true);
439
440 std::cout << "--> setBlock test ended with success." <<std::endl;
441 }
442
testSetBlock2()443 void SimpleMatrixTest::testSetBlock2()
444 {
445 std::cout << "--> Test: testSetBlock2." <<std::endl;
446 // Copy of a Simple into a sub-block of Simple
447 SP::SimpleMatrix MOut(new SimpleMatrix(10, 10));
448
449 SP::SiconosMatrix MIn(new SimpleMatrix(5, 5));
450 for(unsigned int i = 0; i < 5; ++i)
451 for(unsigned int j = 0 ; j < 5; ++j)
452 (*MIn)(i, j) = i + j;
453
454 MOut->setBlock(2, 3, *MIn);
455
456 for(unsigned int i = 2; i < 7; ++i)
457 for(unsigned int j = 3 ; j < 8; ++j)
458 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSetBlock2: ", fabs((*MOut)(i, j) - (*MIn)(i - 2, j - 3)) < tol, true);
459
460 // Copy of a Block into a sub-block of Simple
461
462 MIn.reset(new BlockMatrix(m4, m4, m4, m4));
463 MOut->setBlock(2, 3, *MIn);
464
465 for(unsigned int i = 2; i < 6; ++i)
466 for(unsigned int j = 3 ; j < 7; ++j)
467 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSetBlock2: ", fabs((*MOut)(i, j) - (*MIn)(i - 2, j - 3)) < tol, true);
468
469 std::cout << "--> setBlock2 test ended with success." <<std::endl;
470 }
471
472
testGetSetRowCol()473 void SimpleMatrixTest::testGetSetRowCol()
474 {
475 std::cout << "--> Test: get, set Row and Col." <<std::endl;
476
477 SP::SiconosVector vIn(new SiconosVector(10, 1.2));
478 SP::BlockVector vBIn(new BlockVector());
479 SP::SiconosVector v1(new SiconosVector(3, 2));
480 SP::SiconosVector v2(new SiconosVector(5, 3));
481 SP::SiconosVector v3(new SiconosVector(2, 4));
482 vBIn->insertPtr(v1);
483 vBIn->insertPtr(v2);
484 vBIn->insertPtr(v3);
485
486 // Set row with a SiconosVector
487 C->setRow(4, *vIn);
488 for(unsigned int i = 0; i < C->size(1); ++i)
489 CPPUNIT_ASSERT_EQUAL_MESSAGE("testGetSetRowCol : ", fabs((*C)(4, i) - 1.2) < tol, true);
490
491 // Set col with a SiconosVector
492 C->setCol(4, *vIn);
493 for(unsigned int i = 0; i < C->size(0); ++i)
494 CPPUNIT_ASSERT_EQUAL_MESSAGE("testGetSetRowCol : ", fabs((*C)(i, 4) - 1.2) < tol, true);
495
496 // C->setCol(4, *vBIn);
497 // for (unsigned int i = 0; i< C->size(1); ++i)
498 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testGetSetRowCol : ", fabs((*C)(4,i)- (*vBIn)(i)) < tol, true);
499
500 *C = *A; //reset C
501 vIn->zero();
502 vBIn->zero();
503 // get row and copy it into a SiconosVector
504 C->getRow(4, *vIn);
505 for(unsigned int i = 0; i < C->size(1); ++i)
506 CPPUNIT_ASSERT_EQUAL_MESSAGE("testGetSetRowCol : ", fabs((*C)(4, i) - (*vIn)(i)) < tol, true);
507
508 // get col and copy it into a SiconosVector
509 C->getCol(4, *vIn);
510 for(unsigned int i = 0; i < C->size(0); ++i)
511 CPPUNIT_ASSERT_EQUAL_MESSAGE("testGetSetRowCol : ", fabs((*C)(i, 4) - (*vIn)(i)) < tol, true);
512
513 std::cout << "--> get, set Row and Col tests ended with success." <<std::endl;
514 }
515
testTrans()516 void SimpleMatrixTest::testTrans()
517 {
518 std::cout << "--> Test: trans." <<std::endl;
519
520 // Transpose in place ...
521 SP::SimpleMatrix ref(new SimpleMatrix(*D));
522 SP::SimpleMatrix tRef(new SimpleMatrix(*ref));
523
524 tRef->trans();
525 for(unsigned int i = 0; i < ref->size(0); ++i)
526 for(unsigned int j = 0 ; j < ref->size(1); ++j)
527 if(i == j)
528 CPPUNIT_ASSERT_EQUAL_MESSAGE("testTrans: ", (*tRef)(i, j) == (*ref)(i, j), true);
529 else
530 CPPUNIT_ASSERT_EQUAL_MESSAGE("testTrans: ", (*tRef)(i, j) == (*ref)(j, i), true);
531
532 // Transpose of another matrix ...
533 // Dense
534 tRef->zero();
535
536 tRef->trans(*ref);
537 for(unsigned int i = 0; i < ref->size(0); ++i)
538 for(unsigned int j = 0 ; j < ref->size(1); ++j)
539 if(i == j)
540 CPPUNIT_ASSERT_EQUAL_MESSAGE("testTrans: ", (*tRef)(i, j) == (*ref)(i, j), true);
541 else
542 CPPUNIT_ASSERT_EQUAL_MESSAGE("testTrans: ", (*tRef)(i, j) == (*ref)(j, i), true);
543
544 // Sym
545 ref.reset(new SimpleMatrix(*S));
546 tRef.reset(new SimpleMatrix(*ref));
547 tRef->trans(*ref);
548 CPPUNIT_ASSERT_EQUAL_MESSAGE("testTrans: ", (*tRef) == (*ref), true);
549 // Sparse
550 ref.reset(new SimpleMatrix(*SP));
551 tRef.reset(new SimpleMatrix(*ref));
552 tRef->trans(*ref);
553 // for(unsigned int i = 0; i<ref->size(0); ++i)
554 // {
555 // for(unsigned int j = 0 ; j< ref->size(1); ++j)
556 // if(i==j)
557 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testTrans: ", (*tRef)(i,j) == (*ref)(i,j) , true);
558 // else
559 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testTrans: ", (*tRef)(i,j) == (*ref)(j,i) , true);
560 // }
561 // Banded
562 // ref.reset(new SimpleMatrix(*Band);
563 // tRef.reset(new SimpleMatrix(*ref);
564 // *tRef = trans(*ref);
565 // for(unsigned int i = 0; i<ref->size(0); ++i)
566 // for(unsigned int j = 0 ; j< ref->size(1); ++j)
567 // if(i==j)
568 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testTrans: ", (*tRef)(i,j) == (*ref)(i,j) , true);
569 // else
570 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testTrans: ", (*tRef)(i,j) == (*ref)(j,i) , true);
571 std::cout << "--> test trans ended with success." <<std::endl;
572 }
573
574
testAssignment0()575 void SimpleMatrixTest::testAssignment0()
576 {
577 std::cout << "--> Test: assignment0." <<std::endl;
578
579 // Simple = Simple
580
581 SP::SimpleMatrix ref(new SimpleMatrix(*D));
582 SP::SiconosMatrix tRef(new SimpleMatrix(*SicM));
583 // Dense = any type:
584 *tRef = *ref;
585 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef) == (*ref), true);
586 ref.reset(new SimpleMatrix(*T));
587 SP::SiconosMatrix tRef3(new SimpleMatrix(3, 3));
588 *tRef3 = *ref;
589 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef3) == (*ref), true);
590
591 ref.reset(new SimpleMatrix(*S));
592 *tRef3 = *ref;
593 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef3) == (*ref), true);
594
595 SP::SiconosMatrix tRef4(new SimpleMatrix(4, 4));
596 ref.reset(new SimpleMatrix(*SP));
597 *tRef4 = *ref;
598 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef4) == (*ref), true);
599
600
601 ref.reset(new SimpleMatrix(*Band));
602 *tRef4 = *ref;
603 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef4) == (*ref), true);
604 ref.reset(new SimpleMatrix(*Z));
605 *tRef3 = *ref;
606 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef3) == (*ref), true);
607
608 ref.reset(new SimpleMatrix(*I));
609 *tRef3 = *ref;
610 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef3) == (*ref), true);
611
612 // Triang = Triang, Zero or Identity
613 ref.reset(new SimpleMatrix(*T));
614 tRef.reset(new SimpleMatrix(*T));
615 tRef->zero();
616 *tRef = *ref;
617 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef) == (*ref), true);
618
619 ref.reset(new SimpleMatrix(*Z));
620 *tRef = *ref;
621 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef) == (*ref), true);
622
623 ref.reset(new SimpleMatrix(*I));
624 *tRef = *ref;
625 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef) == (*ref), true);
626 // Sym = Sym, Zero or Id
627 ref.reset(new SimpleMatrix(*S));
628 tRef.reset(new SimpleMatrix(*S));
629 tRef->zero();
630 *tRef = *ref;
631 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef) == (*ref), true);
632 ref.reset(new SimpleMatrix(*Z));
633 *tRef = *ref;
634 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef) == (*ref), true);
635 ref.reset(new SimpleMatrix(*I));
636 *tRef = *ref;
637 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef) == (*ref), true);
638
639 // Sparse => Sparse or Zero
640 ref.reset(new SimpleMatrix(*SP));
641 tRef.reset(new SimpleMatrix(*SP));
642 tRef->zero();
643 *tRef = *ref;
644 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef) == (*ref), true);
645
646 ref.reset(new SimpleMatrix(*Z2));
647 *tRef = *ref;
648
649 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef) == (*ref), true);
650
651 // // Sparse coordinate => Sparse
652 ref.reset(new SimpleMatrix(*SP_coor));
653 //ref->display();
654 tRef.reset(new SimpleMatrix(*SP));
655 tRef->zero();
656 *tRef = *ref;
657 //tRef->displayExpert();
658 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef) == (*ref), true);
659
660
661 // Banded = Banded, Id or Zero
662 ref.reset(new SimpleMatrix(*Band));
663 tRef.reset(new SimpleMatrix(*Band));
664 tRef->zero();
665 *tRef = *ref;
666 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef) == (*ref), true);
667 ref.reset(new SimpleMatrix(*Z2));
668 *tRef = *ref;
669 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef) == (*ref), true);
670 ref.reset(new SimpleMatrix(*I2));
671 *tRef = *ref;
672 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment0: ", (*tRef) == (*ref), true);
673
674
675
676 std::cout << "--> test assignment0 ended with success." <<std::endl;
677 }
678
testAssignment1()679 void SimpleMatrixTest::testAssignment1()
680 {
681 std::cout << "--> Test: assignment1." <<std::endl;
682
683 // Simple = Siconos(Block)
684
685 *C = *Ab;
686 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment1: ", (*C) == (*Ab), true);
687 std::cout << "--> test assignment1 ended with success." <<std::endl;
688 }
689
testAssignment2()690 void SimpleMatrixTest::testAssignment2()
691 {
692 std::cout << "--> Test: assignment2." <<std::endl;
693
694 // Simple = Siconos(Simple)
695
696 SP::SiconosMatrix ref(new SimpleMatrix(*D));
697 SP::SiconosMatrix tRef(new SimpleMatrix(*SicM));
698 // Dense = any type:
699 *tRef = *ref;
700 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef) == (*ref), true);
701
702 ref.reset(new SimpleMatrix(*T));
703 SP::SiconosMatrix tRef3(new SimpleMatrix(3, 3));
704 *tRef3 = *ref;
705 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef3) == (*ref), true);
706
707 ref.reset(new SimpleMatrix(*S));
708 *tRef3 = *ref;
709 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef3) == (*ref), true);
710
711 SP::SiconosMatrix tRef4(new SimpleMatrix(4, 4));
712 ref.reset(new SimpleMatrix(*SP));
713 *tRef4 = *ref;
714 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef4) == (*ref), true);
715
716 ref.reset(new SimpleMatrix(*Band));
717 *tRef4 = *ref;
718 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef4) == (*ref), true);
719 ref.reset(new SimpleMatrix(*Z));
720 *tRef3 = *ref;
721 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef3) == (*ref), true);
722
723 ref.reset(new SimpleMatrix(*I));
724 *tRef3 = *ref;
725 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef3) == (*ref), true);
726 // Triang = Triang, Zero or Identity
727 ref.reset(new SimpleMatrix(*T));
728 tRef.reset(new SimpleMatrix(*T));
729 tRef->zero();
730 *tRef = *ref;
731 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef) == (*ref), true);
732
733 ref.reset(new SimpleMatrix(*Z));
734 *tRef = *ref;
735 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef) == (*ref), true);
736
737 ref.reset(new SimpleMatrix(*I));
738 *tRef = *ref;
739 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef) == (*ref), true);
740 // Sym = Sym, Zero or Id
741 ref.reset(new SimpleMatrix(*S));
742 tRef.reset(new SimpleMatrix(*S));
743 tRef->zero();
744 *tRef = *ref;
745 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef) == (*ref), true);
746 ref.reset(new SimpleMatrix(*Z));
747 *tRef = *ref;
748 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef) == (*ref), true);
749
750 ref.reset(new SimpleMatrix(*I));
751 *tRef = *ref;
752 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef) == (*ref), true);
753 // Sparse = Sparse or Zero
754 ref.reset(new SimpleMatrix(*SP));
755 tRef.reset(new SimpleMatrix(*SP));
756 tRef->zero();
757 *tRef = *ref;
758 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef) == (*ref), true);
759
760 ref.reset(new SimpleMatrix(*Z2));
761 *tRef = *ref;
762 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef) == (*ref), true);
763 // Banded = Banded, Id or Zero
764 ref.reset(new SimpleMatrix(*Band));
765 tRef.reset(new SimpleMatrix(*Band));
766 tRef->zero();
767 *tRef = *ref;
768 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef) == (*ref), true);
769
770 ref.reset(new SimpleMatrix(*Z2));
771 *tRef = *ref;
772 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef) == (*ref), true);
773 ref.reset(new SimpleMatrix(*I2));
774 *tRef = *ref;
775 CPPUNIT_ASSERT_EQUAL_MESSAGE("testAssignment2: ", (*tRef) == (*ref), true);
776 std::cout << "--> test assignment2 ended with success." <<std::endl;
777 }
778
testOperators1()779 void SimpleMatrixTest::testOperators1()
780 {
781 std::cout << "--> Test: operators1." <<std::endl;
782 //+=, -=, *=, /=
783
784 SP::SiconosMatrix tmp(new SimpleMatrix(*D));
785 // Dense *=, /=
786 double a = 2.2;
787 int a1 = 2;
788 *tmp *= a;
789 for(unsigned int i = 0; i < tmp->size(0); ++i)
790 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
791 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*tmp)(i, j) - a * (*D)(i, j)) < tol, true);
792
793 *tmp *= a1;
794 for(unsigned int i = 0; i < tmp->size(0); ++i)
795 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
796 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*tmp)(i, j) - a * a1 * (*D)(i, j)) < tol, true);
797
798 *tmp /= a;
799 for(unsigned int i = 0; i < tmp->size(0); ++i)
800 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
801 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*tmp)(i, j) - a1 * (*D)(i, j)) < tol, true);
802
803 *tmp /= a1;
804 for(unsigned int i = 0; i < tmp->size(0); ++i)
805 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
806 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*tmp)(i, j) - (*D)(i, j)) < tol, true);
807
808 // Dense +=, -= Dense
809
810 *tmp += *SicM;
811 for(unsigned int i = 0; i < tmp->size(0); ++i)
812 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
813 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*tmp)(i, j) - (*SicM)(i, j) - (*D)(i, j)) < tol, true);
814
815 *tmp -= *SicM;
816 for(unsigned int i = 0; i < tmp->size(0); ++i)
817 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
818 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*tmp)(i, j) - (*D)(i, j)) < tol, true);
819
820 // Dense +=, -= Block
821 C->zero();
822 *C += *Ab;
823 *C += *Ab;
824 for(unsigned int i = 0; i < C->size(0); ++i)
825 for(unsigned int j = 0 ; j < C->size(1); ++j)
826 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*C)(i, j) - 2 * (*Ab)(i, j)) < tol, true);
827 *C -= *Ab;
828 for(unsigned int i = 0; i < C->size(0); ++i)
829 for(unsigned int j = 0 ; j < C->size(1); ++j)
830 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*C)(i, j) - (*Ab)(i, j)) < tol, true);
831
832 std::cout << "--> test operators1 ended with success." <<std::endl;
833 }
834
testOperators2()835 void SimpleMatrixTest::testOperators2()
836 {
837 std::cout << "--> Test: operators2." <<std::endl;
838 // +=, -=, *=, /= triangular
839 SP::SiconosMatrix tmp(new SimpleMatrix(*T));
840 SP::SiconosMatrix tmp2(new SimpleMatrix(*T));
841 *tmp += *tmp2;
842 for(unsigned int i = 0; i < tmp->size(0); ++i)
843 for(unsigned int j = i ; j < tmp->size(1); ++j)
844 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", (*tmp)(i, j) == 2.0 * (*T)(i, j), true);
845
846 int mult = 2;
847 double mult0 = 2.2;
848 *tmp *= mult0;
849 for(unsigned int i = 0; i < tmp->size(0); ++i)
850 for(unsigned int j = i ; j < tmp->size(1); ++j)
851 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult0 * (*T)(i, j)) < tol, true);
852
853 *tmp *= mult;
854 for(unsigned int i = 0; i < tmp->size(0); ++i)
855 for(unsigned int j = i ; j < tmp->size(1); ++j)
856 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult * mult0 * (*T)(i, j)) < tol, true);
857
858 *tmp /= mult;
859 for(unsigned int i = 0; i < tmp->size(0); ++i)
860 for(unsigned int j = i ; j < tmp->size(1); ++j)
861 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult0 * (*T)(i, j)) < tol, true);
862
863 *tmp /= mult0;
864 for(unsigned int i = 0; i < tmp->size(0); ++i)
865 for(unsigned int j = i ; j < tmp->size(1); ++j)
866 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", (*tmp)(i, j) == 2 * (*T)(i, j), true);
867
868 *tmp -= *tmp2;
869 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(tmp->getTriang() - *T) == 0, true);
870 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", tmp->num() == Siconos::TRIANGULAR, true);
871
872 std::cout << "--> test operators2 ended with success." <<std::endl;
873 }
874
testOperators3()875 void SimpleMatrixTest::testOperators3()
876 {
877 std::cout << "--> Test: operators3." <<std::endl;
878 // +=, -=, *=, /= Symmetric
879 SP::SiconosMatrix tmp(new SimpleMatrix(*S));
880 SP::SiconosMatrix tmp2(new SimpleMatrix(*S));
881 *tmp += *tmp2;
882 for(unsigned int i = 0; i < tmp->size(0); ++i)
883 for(unsigned int j = i ; j < tmp->size(1); ++j)
884 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", (*tmp)(i, j) == 2.0 * (*S)(i, j), true);
885
886 int mult = 2;
887 double mult0 = 2.2;
888 *tmp *= mult0;
889 for(unsigned int i = 0; i < tmp->size(0); ++i)
890 for(unsigned int j = i ; j < tmp->size(1); ++j)
891 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult0 * (*S)(i, j)) < tol, true);
892
893 *tmp *= mult;
894 for(unsigned int i = 0; i < tmp->size(0); ++i)
895 for(unsigned int j = i ; j < tmp->size(1); ++j)
896 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult * mult0 * (*S)(i, j)) < tol, true);
897
898 *tmp /= mult;
899 for(unsigned int i = 0; i < tmp->size(0); ++i)
900 for(unsigned int j = i ; j < tmp->size(1); ++j)
901 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult0 * (*S)(i, j)) < tol, true);
902
903 *tmp /= mult0;
904 for(unsigned int i = 0; i < tmp->size(0); ++i)
905 for(unsigned int j = i ; j < tmp->size(1); ++j)
906 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", (*tmp)(i, j) == 2 * (*S)(i, j), true);
907
908 *tmp -= *tmp2;
909 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(tmp->getSym() - *S) == 0, true);
910 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", tmp->num() == Siconos::SYMMETRIC, true);
911
912 std::cout << "--> test operators3 ended with success." <<std::endl;
913 }
914
testOperators4()915 void SimpleMatrixTest::testOperators4()
916 {
917 std::cout << "--> Test: operators4." <<std::endl;
918 // +=, -=, *=, /= sparse
919 SP::SiconosMatrix tmp(new SimpleMatrix(*SP));
920 SP::SiconosMatrix tmp2(new SimpleMatrix(*SP));
921 SP::SiconosMatrix tmp3(new SimpleMatrix(*T2));
922
923 SP::SiconosMatrix tmp4(new SimpleMatrix(*Band));
924 SP::SiconosMatrix tmp5(new SimpleMatrix(*S2));
925
926 *tmp += *tmp2;
927 for(unsigned int i = 0; i < tmp->size(0); ++i)
928 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
929 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * (*SP)(i, j)) < tol, true);
930
931 int mult = 2;
932 double mult0 = 2.2;
933 *tmp *= mult0;
934 for(unsigned int i = 0; i < tmp->size(0); ++i)
935 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
936 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult0 * (*SP)(i, j)) < tol, true);
937
938 *tmp *= mult;
939 for(unsigned int i = 0; i < tmp->size(0); ++i)
940 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
941 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult * mult0 * (*SP)(i, j)) < tol, true);
942
943 *tmp /= mult;
944 for(unsigned int i = 0; i < tmp->size(0); ++i)
945 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
946 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult0 * (*SP)(i, j)) < tol, true);
947
948 *tmp /= mult0;
949 for(unsigned int i = 0; i < tmp->size(0); ++i)
950 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
951 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2 * (*SP)(i, j)) < tol, true);
952
953 *tmp -= *tmp2;
954 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(tmp->getSparse() - *SP) == 0, true);
955 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", tmp->num() == Siconos::SPARSE, true);
956
957 // += -= a triangular
958 *tmp += *tmp3;
959 for(unsigned int i = 0; i < tmp->size(0); ++i)
960 {
961 for(unsigned int j = 0; j < i; ++j)
962 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP)(i, j)) < tol, true);
963 for(unsigned int j = i ; j < tmp->size(1); ++j)
964 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP)(i, j) - (*tmp3)(i, j)) < tol, true);
965 }
966
967 *tmp -= *tmp3;
968 for(unsigned int i = 0; i < tmp->size(0); ++i)
969 for(unsigned int j = 0; j < tmp->size(1); ++j)
970 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP)(i, j)) < tol, true);
971
972 // += -= a banded
973 *tmp -= *tmp;
974 *tmp += *tmp4;
975 for(signed i = 0; i < signed(Band->size1()); ++ i)
976 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
977 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*Band)(i, j)) < tol, true);
978
979 *tmp -= *tmp4;
980 for(unsigned int i = 0; i < tmp->size(0); ++i)
981 for(unsigned int j = 0; j < tmp->size(1); ++j)
982 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", (*tmp)(i, j) < tol, true);
983
984 // += -= a sym
985
986 *tmp += *tmp5;
987 for(unsigned int i = 0; i < tmp->size(0); ++i)
988 {
989 for(unsigned int j = 0; j < i; ++j)
990 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP)(i, j) - (*tmp5)(j, i)) < tol, true);
991 for(unsigned int j = i ; j < tmp->size(1); ++j)
992 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP)(i, j) - (*tmp5)(i, j)) < tol, true);
993 }
994
995 *tmp -= *tmp5;
996 for(unsigned int i = 0; i < tmp->size(0); ++i)
997 for(unsigned int j = 0; j < tmp->size(1); ++j)
998 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP)(i, j)) < tol, true);
999
1000 std::cout << "--> test operators4 ended with success." <<std::endl;
1001 }
testOperators4bis()1002 void SimpleMatrixTest::testOperators4bis()
1003 {
1004 std::cout << "--> Test: operators4bis." <<std::endl;
1005 // +=, -=, *=, /= sparse
1006 SP::SiconosMatrix tmp(new SimpleMatrix(*SP_coor));
1007 SP::SiconosMatrix tmp2(new SimpleMatrix(*SP_coor));
1008 SP::SiconosMatrix tmp3(new SimpleMatrix(*T2));
1009
1010 SP::SiconosMatrix tmp4(new SimpleMatrix(*Band));
1011 SP::SiconosMatrix tmp5(new SimpleMatrix(*S2));
1012
1013 SP::SiconosMatrix tmp6(new SimpleMatrix(*SP));
1014
1015
1016 *tmp += *tmp2;
1017 for(unsigned int i = 0; i < tmp->size(0); ++i)
1018 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
1019 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * (*SP_coor)(i, j)) < tol, true);
1020
1021 int mult = 2;
1022 double mult0 = 2.2;
1023 *tmp *= mult0;
1024 for(unsigned int i = 0; i < tmp->size(0); ++i)
1025 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
1026 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult0 * (*SP_coor)(i, j)) < tol, true);
1027
1028 *tmp *= mult;
1029 for(unsigned int i = 0; i < tmp->size(0); ++i)
1030 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
1031 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult * mult0 * (*SP_coor)(i, j)) < tol, true);
1032
1033 *tmp /= mult;
1034 for(unsigned int i = 0; i < tmp->size(0); ++i)
1035 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
1036 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult0 * (*SP_coor)(i, j)) < tol, true);
1037
1038 *tmp /= mult0;
1039 for(unsigned int i = 0; i < tmp->size(0); ++i)
1040 for(unsigned int j = 0 ; j < tmp->size(1); ++j)
1041 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2 * (*SP_coor)(i, j)) < tol, true);
1042
1043 *tmp -= *tmp2;
1044 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(tmp->getSparseCoordinate() - *SP_coor) == 0, true);
1045 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", tmp->num() == SPARSE_COORDINATE, true);
1046
1047 // += -= a triangular
1048 *tmp += *tmp3;
1049 for(unsigned int i = 0; i < tmp->size(0); ++i)
1050 {
1051 for(unsigned int j = 0; j < i; ++j)
1052 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP_coor)(i, j)) < tol, true);
1053 for(unsigned int j = i ; j < tmp->size(1); ++j)
1054 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP_coor)(i, j) - (*tmp3)(i, j)) < tol, true);
1055 }
1056
1057 *tmp -= *tmp3;
1058 for(unsigned int i = 0; i < tmp->size(0); ++i)
1059 for(unsigned int j = 0; j < tmp->size(1); ++j)
1060 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP_coor)(i, j)) < tol, true);
1061
1062 // += -= a banded
1063 *tmp -= *tmp;
1064 *tmp += *tmp4;
1065 for(signed i = 0; i < signed(Band->size1()); ++ i)
1066 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1067 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*Band)(i, j)) < tol, true);
1068
1069 *tmp -= *tmp4;
1070 for(unsigned int i = 0; i < tmp->size(0); ++i)
1071 for(unsigned int j = 0; j < tmp->size(1); ++j)
1072 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", (*tmp)(i, j) < tol, true);
1073
1074 // += -= a sym
1075
1076 *tmp += *tmp5;
1077 for(unsigned int i = 0; i < tmp->size(0); ++i)
1078 {
1079 for(unsigned int j = 0; j < i; ++j)
1080 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP_coor)(i, j) - (*tmp5)(j, i)) < tol, true);
1081 for(unsigned int j = i ; j < tmp->size(1); ++j)
1082 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP_coor)(i, j) - (*tmp5)(i, j)) < tol, true);
1083 }
1084
1085 *tmp -= *tmp5;
1086 for(unsigned int i = 0; i < tmp->size(0); ++i)
1087 for(unsigned int j = 0; j < tmp->size(1); ++j)
1088 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP_coor)(i, j)) < tol, true);
1089
1090 // += -= a sparse
1091
1092 *tmp += *tmp6;
1093 for(unsigned int i = 0; i < tmp->size(0); ++i)
1094 {
1095 for(unsigned int j = 0; j < i; ++j)
1096 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP_coor)(i, j) - (*tmp6)(j, i)) < tol, true);
1097 for(unsigned int j = i ; j < tmp->size(1); ++j)
1098 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP_coor)(i, j) - (*tmp6)(i, j)) < tol, true);
1099 }
1100
1101 *tmp -= *tmp6;
1102 for(unsigned int i = 0; i < tmp->size(0); ++i)
1103 for(unsigned int j = 0; j < tmp->size(1); ++j)
1104 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - (*SP_coor)(i, j)) < tol, true);
1105
1106 std::cout << "--> test operators4bis ended with success." <<std::endl;
1107 }
testOperators5()1108 void SimpleMatrixTest::testOperators5()
1109 {
1110 std::cout << "--> Test: operators5." <<std::endl;
1111 // +=, -=, *=, /= banded
1112 SP::SiconosMatrix tmp(new SimpleMatrix(*Band));
1113 SP::SiconosMatrix tmp2(new SimpleMatrix(*Band));
1114 *tmp += *tmp2;
1115 for(signed i = 0; i < signed(Band->size1()); ++ i)
1116 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1117 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", (*tmp)(i, j) == 2.0 * (*Band)(i, j), true);
1118
1119 int mult = 2;
1120 double mult0 = 2.2;
1121 *tmp *= mult0;
1122 for(signed i = 0; i < signed(Band->size1()); ++ i)
1123 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1124 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult0 * (*Band)(i, j)) < tol, true);
1125
1126 *tmp *= mult;
1127 for(signed i = 0; i < signed(Band->size1()); ++ i)
1128 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1129 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult * mult0 * (*Band)(i, j)) < tol, true);
1130
1131 *tmp /= mult;
1132 for(signed i = 0; i < signed(Band->size1()); ++ i)
1133 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1134 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", ((*tmp)(i, j) - 2.0 * mult0 * (*Band)(i, j)) < tol, true);
1135
1136 *tmp /= mult0;
1137 for(signed i = 0; i < signed(Band->size1()); ++ i)
1138 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1139 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", (*tmp)(i, j) == 2 * (*Band)(i, j), true);
1140
1141 *tmp -= *tmp2;
1142 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(tmp->getBanded() - *Band) == 0, true);
1143 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", tmp->num() == Siconos::BANDED, true);
1144
1145 std::cout << "--> test operators5 ended with success." <<std::endl;
1146 }
1147
testOperators6()1148 void SimpleMatrixTest::testOperators6()
1149 {
1150 std::cout << "--> Test: operator6." <<std::endl;
1151
1152 // ============= C = A + B =============
1153
1154 // Dense = Dense + Dense
1155 C->zero();
1156 *C = *A + *B;
1157 for(unsigned int i = 0; i < C->size(0); ++i)
1158 for(unsigned int j = 0 ; j < C->size(1); ++j)
1159 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*C)(i, j) - (*A)(i, j) - (*B)(i, j)) < tol, true);
1160
1161 C->zero();
1162 // Dense = Dense + Block
1163 *C = *A + *Ab;
1164 for(unsigned int i = 0; i < C->size(0); ++i)
1165 for(unsigned int j = 0 ; j < C->size(1); ++j)
1166 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*C)(i, j) - (*A)(i, j) - (*Ab)(i, j)) < tol, true);
1167
1168 C->zero();
1169 // Dense = Block + Dense
1170 *C = *Ab + *A;
1171 for(unsigned int i = 0; i < C->size(0); ++i)
1172 for(unsigned int j = 0 ; j < C->size(1); ++j)
1173 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*C)(i, j) - (*A)(i, j) - (*Ab)(i, j)) < tol, true);
1174
1175 C->zero();
1176
1177 // Dense = Block + Block
1178 *C = *Ab + *Ab;
1179 for(unsigned int i = 0; i < C->size(0); ++i)
1180 for(unsigned int j = 0 ; j < C->size(1); ++j)
1181 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*C)(i, j) - (*Ab)(i, j) - (*Ab)(i, j)) < tol, true);
1182
1183 C->zero();
1184
1185 // Block = Dense + Dense
1186 Cb->zero();
1187 *Cb = *A + *B;
1188 for(unsigned int i = 0; i < Cb->size(0); ++i)
1189 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1190 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*Cb)(i, j) - (*A)(i, j) - (*B)(i, j)) < tol, true);
1191
1192 // Block = Dense + Block
1193
1194 Cb->zero();
1195 *Cb = *A + *Ab;
1196 for(unsigned int i = 0; i < Cb->size(0); ++i)
1197 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1198 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*Cb)(i, j) - (*A)(i, j) - (*Ab)(i, j)) < tol, true);
1199
1200
1201 // Block = Block + Dense
1202
1203 Cb->zero();
1204 *Cb = *Ab + *A;
1205 for(unsigned int i = 0; i < Cb->size(0); ++i)
1206 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1207 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*Cb)(i, j) - (*Ab)(i, j) - (*A)(i, j)) < tol, true);
1208
1209
1210 // Block = Block + Block
1211
1212 Cb->zero();
1213 *Cb = *Ab + *Bb;
1214 for(unsigned int i = 0; i < Cb->size(0); ++i)
1215 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1216 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*Cb)(i, j) - (*Ab)(i, j) - (*Bb)(i, j)) < tol, true);
1217 Cb->zero();
1218
1219 // ============= C = A - B =============
1220
1221 // Dense = Dense - Dense
1222 C->zero();
1223 *C = *A - *B;
1224 for(unsigned int i = 0; i < C->size(0); ++i)
1225 for(unsigned int j = 0 ; j < C->size(1); ++j)
1226 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*C)(i, j) - (*A)(i, j) + (*B)(i, j)) < tol, true);
1227
1228 C->zero();
1229 // Dense = Dense - Block
1230 *C = *A - *Ab;
1231 for(unsigned int i = 0; i < C->size(0); ++i)
1232 for(unsigned int j = 0 ; j < C->size(1); ++j)
1233 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*C)(i, j) - (*A)(i, j) + (*Ab)(i, j)) < tol, true);
1234
1235 C->zero();
1236 // Dense = Block - Dense
1237 *C = *Ab - *A;
1238 for(unsigned int i = 0; i < C->size(0); ++i)
1239 for(unsigned int j = 0 ; j < C->size(1); ++j)
1240 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*C)(i, j) + (*A)(i, j) - (*Ab)(i, j)) < tol, true);
1241
1242 C->zero();
1243
1244 // Dense = Block - Block
1245 *C = *Ab - *Bb;
1246 for(unsigned int i = 0; i < C->size(0); ++i)
1247 for(unsigned int j = 0 ; j < C->size(1); ++j)
1248 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*C)(i, j) - (*Ab)(i, j) + (*Bb)(i, j)) < tol, true);
1249
1250 C->zero();
1251
1252 // Block = Dense - Dense
1253 Cb->zero();
1254 *Cb = *A - *B;
1255 for(unsigned int i = 0; i < Cb->size(0); ++i)
1256 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1257 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*Cb)(i, j) - (*A)(i, j) + (*B)(i, j)) < tol, true);
1258
1259 // Block = Dense - Block
1260
1261 Cb->zero();
1262 *Cb = *A - *Ab;
1263 for(unsigned int i = 0; i < Cb->size(0); ++i)
1264 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1265 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*Cb)(i, j) - (*A)(i, j) + (*Ab)(i, j)) < tol, true);
1266
1267
1268 // Block = Block - Dense
1269
1270 Cb->zero();
1271 *Cb = *Ab - *A;
1272 for(unsigned int i = 0; i < Cb->size(0); ++i)
1273 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1274 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*Cb)(i, j) - (*Ab)(i, j) + (*A)(i, j)) < tol, true);
1275
1276
1277 // Block = Block - Block
1278
1279 Cb->zero();
1280 *Cb = *Ab - *Bb;
1281 for(unsigned int i = 0; i < Cb->size(0); ++i)
1282 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1283 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6: ", fabs((*Cb)(i, j) - (*Ab)(i, j) + (*Bb)(i, j)) < tol, true);
1284 Cb->zero();
1285 std::cout << "--> test operators6 ended with success." <<std::endl;
1286 }
1287
testOperators6Bis()1288 void SimpleMatrixTest::testOperators6Bis()
1289 {
1290 std::cout << "--> Test: operator6Bis." <<std::endl;
1291
1292 // ============= C = A + B =============
1293
1294 // Dense = Dense + Dense
1295 C->zero();
1296 add(*A, *B, *C);
1297 for(unsigned int i = 0; i < C->size(0); ++i)
1298 for(unsigned int j = 0 ; j < C->size(1); ++j)
1299 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*C)(i, j) - (*A)(i, j) - (*B)(i, j)) < tol, true);
1300
1301 C->zero();
1302 // Dense = Dense + Block
1303 add(*A, *Ab, *C);
1304 for(unsigned int i = 0; i < C->size(0); ++i)
1305 for(unsigned int j = 0 ; j < C->size(1); ++j)
1306 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*C)(i, j) - (*A)(i, j) - (*Ab)(i, j)) < tol, true);
1307
1308 C->zero();
1309 // Dense = Block + Dense
1310 add(*Ab, *A, *C);
1311 for(unsigned int i = 0; i < C->size(0); ++i)
1312 for(unsigned int j = 0 ; j < C->size(1); ++j)
1313 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*C)(i, j) - (*A)(i, j) - (*Ab)(i, j)) < tol, true);
1314
1315 C->zero();
1316
1317 // Dense = Block + Block
1318 add(*Ab, *Ab, *C);
1319 for(unsigned int i = 0; i < C->size(0); ++i)
1320 for(unsigned int j = 0 ; j < C->size(1); ++j)
1321 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*C)(i, j) - (*Ab)(i, j) - (*Ab)(i, j)) < tol, true);
1322
1323 C->zero();
1324
1325 // Block = Dense + Dense
1326 Cb->zero();
1327 add(*A, *B, *Cb);
1328 for(unsigned int i = 0; i < Cb->size(0); ++i)
1329 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1330 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*Cb)(i, j) - (*A)(i, j) - (*B)(i, j)) < tol, true);
1331
1332 // Block = Dense + Block
1333
1334 Cb->zero();
1335 add(*A, *Ab, *Cb);
1336 for(unsigned int i = 0; i < Cb->size(0); ++i)
1337 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1338 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*Cb)(i, j) - (*A)(i, j) - (*Ab)(i, j)) < tol, true);
1339
1340
1341 // Block = Block + Dense
1342
1343 Cb->zero();
1344 add(*Ab, *A, *Cb);
1345 for(unsigned int i = 0; i < Cb->size(0); ++i)
1346 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1347 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*Cb)(i, j) - (*Ab)(i, j) - (*A)(i, j)) < tol, true);
1348
1349
1350 // Block = Block + Block
1351
1352 Cb->zero();
1353 add(*Ab, *Bb, *Cb);
1354 for(unsigned int i = 0; i < Cb->size(0); ++i)
1355 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1356 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*Cb)(i, j) - (*Ab)(i, j) - (*Bb)(i, j)) < tol, true);
1357 Cb->zero();
1358
1359 // ============= C = A - B =============
1360
1361 // Dense = Dense - Dense
1362 C->zero();
1363 sub(*A, *B, *C);
1364 for(unsigned int i = 0; i < C->size(0); ++i)
1365 for(unsigned int j = 0 ; j < C->size(1); ++j)
1366 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*C)(i, j) - (*A)(i, j) + (*B)(i, j)) < tol, true);
1367
1368 C->zero();
1369 // Dense = Dense - Block
1370 sub(*A, *Ab, *C);
1371 for(unsigned int i = 0; i < C->size(0); ++i)
1372 for(unsigned int j = 0 ; j < C->size(1); ++j)
1373 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*C)(i, j) - (*A)(i, j) + (*Ab)(i, j)) < tol, true);
1374
1375 C->zero();
1376 // Dense = Block - Dense
1377 sub(*Ab, *A, *C);
1378 for(unsigned int i = 0; i < C->size(0); ++i)
1379 for(unsigned int j = 0 ; j < C->size(1); ++j)
1380 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*C)(i, j) + (*A)(i, j) - (*Ab)(i, j)) < tol, true);
1381
1382 C->zero();
1383
1384 // Dense = Block - Block
1385 sub(*Ab, *Bb, *C);
1386 for(unsigned int i = 0; i < C->size(0); ++i)
1387 for(unsigned int j = 0 ; j < C->size(1); ++j)
1388 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*C)(i, j) - (*Ab)(i, j) + (*Bb)(i, j)) < tol, true);
1389
1390 C->zero();
1391
1392 // Block = Dense - Dense
1393 Cb->zero();
1394 sub(*A, *B, *Cb);
1395 for(unsigned int i = 0; i < Cb->size(0); ++i)
1396 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1397 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*Cb)(i, j) - (*A)(i, j) + (*B)(i, j)) < tol, true);
1398
1399 // Block = Dense - Block
1400
1401 Cb->zero();
1402 sub(*A, *Ab, *Cb);
1403 for(unsigned int i = 0; i < Cb->size(0); ++i)
1404 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1405 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*Cb)(i, j) - (*A)(i, j) + (*Ab)(i, j)) < tol, true);
1406
1407
1408 // Block = Block - Dense
1409
1410 Cb->zero();
1411 sub(*Ab, *A, *Cb);
1412 for(unsigned int i = 0; i < Cb->size(0); ++i)
1413 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1414 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*Cb)(i, j) - (*Ab)(i, j) + (*A)(i, j)) < tol, true);
1415
1416
1417 // Block = Block - Block
1418
1419 Cb->zero();
1420 sub(*Ab, *Bb, *Cb);
1421 for(unsigned int i = 0; i < Cb->size(0); ++i)
1422 for(unsigned int j = 0 ; j < Cb->size(1); ++j)
1423 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Bis: ", fabs((*Cb)(i, j) - (*Ab)(i, j) + (*Bb)(i, j)) < tol, true);
1424 Cb->zero();
1425 std::cout << "--> test operators6Bis ended with success." <<std::endl;
1426 }
1427
testOperators6Ter()1428 void SimpleMatrixTest::testOperators6Ter()
1429 {
1430 std::cout << "--> Test: operator6Ter." <<std::endl;
1431
1432 // +, - for non-dense matrices.
1433
1434 // Triang +,-,* Triang
1435 SP::SiconosMatrix tmp(new SimpleMatrix(*T));
1436 SP::SiconosMatrix tmp2(new SimpleMatrix(*T));
1437 SP::SiconosMatrix res(new SimpleMatrix(3, 3, TRIANGULAR));
1438 *res = *tmp + *tmp2;
1439 for(unsigned int i = 0; i < res->size(0); ++i)
1440 for(unsigned int j = i ; j < res->size(1); ++j)
1441 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Ter: ", (*res)(i, j) == ((*T)(i, j) + (*T)(i, j)), true);
1442
1443 *res = *tmp - *tmp2;
1444 for(unsigned int i = 0; i < res->size(0); ++i)
1445 for(unsigned int j = i ; j < res->size(1); ++j)
1446 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Ter: ", (*res)(i, j) == 0, true);
1447
1448 // prod(*tmp, *tmp2, *res);
1449 // prod(*T,*T, *tmp, true);
1450 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Ter: ", norm_inf(res->getTriang() - *tmp) == 0, true);
1451
1452
1453 // Sym +,-,* Sym
1454 tmp.reset(new SimpleMatrix(*S));
1455 tmp2.reset(new SimpleMatrix(*S));
1456 res.reset(new SimpleMatrix(3, 3, SYMMETRIC));
1457 *res = *tmp + *tmp2;
1458 for(unsigned int i = 0; i < res->size(0); ++i)
1459 for(unsigned int j = i ; j < res->size(1); ++j)
1460 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Ter: ", (*res)(i, j) == ((*S)(i, j) + (*S)(i, j)), true);
1461
1462 *res = *tmp - *tmp2;
1463 for(unsigned int i = 0; i < res->size(0); ++i)
1464 for(unsigned int j = i ; j < res->size(1); ++j)
1465 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Ter: ", (*res)(i, j) == 0, true);
1466
1467 // prod(*tmp , *tmp2, *res, true);
1468
1469 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Ter: ", norm_inf(res->getSym() - prod(*S, *S)) == 0, true);
1470
1471
1472 // Sparse +,-,* Sparse
1473 tmp.reset(new SimpleMatrix(*SP));
1474 tmp2.reset(new SimpleMatrix(*SP));
1475 res.reset(new SimpleMatrix(4, 4, Siconos::SPARSE));
1476 *res = *tmp + *tmp2;
1477 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Ter: ", (*res) == (2.0 * (*tmp)), true);
1478
1479 // *res = prod(*tmp , *tmp2);
1480 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Ter: ", norm_inf(*res->sparse() - prod(*SP, *SP)) < tol, true);
1481
1482 *res = *tmp - *tmp2;
1483 tmp->zero();
1484 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Ter: ", (*res) == *tmp, true);
1485
1486 // SparseCoordinate +,-,* SparseCoordinate
1487
1488 tmp.reset(new SimpleMatrix(*SP_coor));
1489 tmp2.reset(new SimpleMatrix(*SP_coor));
1490 res.reset(new SimpleMatrix(4, 4, Siconos::SPARSE_COORDINATE));
1491 *res = *tmp + *tmp2;
1492
1493 // res->displayExpert();
1494 // tmp->displayExpert();
1495 // tmp2->displayExpert();
1496
1497
1498 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Ter: ", (*res) == (2.0 * (*tmp)), true);
1499
1500 // *res = prod(*tmp , *tmp2);
1501 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Ter: ", norm_inf(*res->sparseCoordinate() - prod(*SP_coor, *SP_coor)) < tol, true);
1502
1503 *res = *tmp - *tmp2;
1504 tmp->zero();
1505 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Ter: ", (*res) == *tmp, true);
1506
1507
1508 // Banded +,- Banded
1509 tmp.reset(new SimpleMatrix(*Band));
1510 tmp2.reset(new SimpleMatrix(*Band));
1511 res.reset(new SimpleMatrix(4, 4, BANDED));
1512 *res = *tmp + *tmp2;
1513 for(signed i = 0; i < signed(Band->size1()); ++ i)
1514 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1515 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Ter: ", (*res)(i, j) == ((*Band)(i, j) + (*Band)(i, j)), true);
1516 *res = *tmp - *tmp2;
1517 for(signed i = 0; i < signed(Band->size1()); ++ i)
1518 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1519 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators6Ter: ", (*res)(i, j) == 0, true);
1520
1521 std::cout << "--> test operators6Ter6 ended with success." <<std::endl;
1522 }
1523
testOperators7()1524 void SimpleMatrixTest::testOperators7()
1525 {
1526 std::cout << "--> Test: operator7." <<std::endl;
1527 SP::SiconosMatrix tmp1(new SimpleMatrix(*D));
1528 tmp1->resize(4, 4);
1529 SP::SiconosMatrix tmp2(new SimpleMatrix(*T2));
1530 SP::SiconosMatrix tmp3(new SimpleMatrix(*S2));
1531 SP::SiconosMatrix tmp4(new SimpleMatrix(*SP));
1532 SP::SiconosMatrix tmp5(new SimpleMatrix(*Band));
1533 SP::SiconosMatrix tmp6(new SimpleMatrix(*Z2));
1534 SP::SiconosMatrix tmp7(new SimpleMatrix(*I2));
1535
1536 SP::SiconosMatrix res(new SimpleMatrix(4, 4));
1537
1538 // dense + ...
1539 // ... triang
1540 add(*tmp1, * tmp2, *res);
1541 for(unsigned int i = 0; i < res->size(0); ++i)
1542 {
1543 for(unsigned int j = i ; j < res->size(1); ++j)
1544 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) - (*tmp2)(i, j)) < tol, true);
1545 for(unsigned int j = 0 ; j < i; ++j)
1546 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j)) < tol, true);
1547 }
1548 // ... Sym
1549 add(*tmp1, * tmp3, *res);
1550 for(unsigned int i = 0; i < res->size(0); ++i)
1551 {
1552 for(unsigned int j = i ; j < res->size(1); ++j)
1553 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) - (*tmp3)(i, j)) < tol, true);
1554 for(unsigned int j = 0 ; j < i; ++j)
1555 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) - (*tmp3)(j, i)) < tol, true);
1556 }
1557 // ... Sparse
1558 add(*tmp1, * tmp4, *res);
1559 for(unsigned int i = 0; i < res->size(0); ++i)
1560 for(unsigned int j = 0 ; j < res->size(1); ++j)
1561 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) - (*tmp4)(i, j)) < tol, true);
1562 // ... Banded
1563 add(*tmp1, * tmp5, *res);
1564 for(signed i = 0; i < signed(Band->size1()); ++ i)
1565 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1566 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) - (*Band)(i, j)) < tol, true);
1567 // Zero
1568 add(*tmp1, * tmp6, *res);
1569 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", *res == *tmp1, true);
1570
1571 // Id
1572 add(*tmp1, * tmp7, *res);
1573 for(unsigned int i = 0; i < res->size(0); ++i)
1574 {
1575 for(unsigned int j = 0 ; j < res->size(1); ++j)
1576 {
1577 if(i == j)
1578 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) - 1) < tol, true);
1579 else
1580 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j)) < tol, true);
1581 }
1582 }
1583
1584 // dense - ...
1585 // ... triangular
1586 sub(*tmp1, * tmp2, *res);
1587 for(unsigned int i = 0; i < res->size(0); ++i)
1588 {
1589 for(unsigned int j = i ; j < res->size(1); ++j)
1590 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) + (*tmp2)(i, j)) < tol, true);
1591 for(unsigned int j = 0 ; j < i; ++j)
1592 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j)) < tol, true);
1593 }
1594 // ... Sym
1595 sub(*tmp1, * tmp3, *res);
1596 for(unsigned int i = 0; i < res->size(0); ++i)
1597 {
1598 for(unsigned int j = i ; j < res->size(1); ++j)
1599 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) + (*tmp3)(i, j)) < tol, true);
1600 for(unsigned int j = 0 ; j < i; ++j)
1601 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) + (*tmp3)(j, i)) < tol, true);
1602 }
1603 // ... Sparse
1604 sub(*tmp1, * tmp4, *res);
1605 for(unsigned int i = 0; i < res->size(0); ++i)
1606 for(unsigned int j = 0 ; j < res->size(1); ++j)
1607 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) + (*tmp4)(i, j)) < tol, true);
1608 // ... Banded
1609 sub(*tmp1, * tmp5, *res);
1610 for(signed i = 0; i < signed(Band->size1()); ++ i)
1611 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1612 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) + (*Band)(i, j)) < tol, true);
1613
1614 // Zero
1615 sub(*tmp1, * tmp6, *res);
1616 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", *res == *tmp1, true);
1617
1618 // Id
1619 sub(*tmp1, * tmp7, *res);
1620 for(unsigned int i = 0; i < res->size(0); ++i)
1621 {
1622 for(unsigned int j = 0 ; j < res->size(1); ++j)
1623 {
1624 if(i == j)
1625 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) + 1) < tol, true);
1626 else
1627 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j)) < tol, true);
1628 }
1629 }
1630 // triang + ...
1631 // ... dense
1632 add(*tmp2, * tmp1, *res);
1633 for(unsigned int i = 0; i < res->size(0); ++i)
1634 {
1635 for(unsigned int j = i ; j < res->size(1); ++j)
1636 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) - (*tmp2)(i, j)) < tol, true);
1637 for(unsigned int j = 0 ; j < i; ++j)
1638 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j)) < tol, true);
1639 }
1640 // ... Sym
1641 add(*tmp2, * tmp3, *res);
1642 for(unsigned int i = 0; i < res->size(0); ++i)
1643 {
1644 for(unsigned int j = i ; j < res->size(1); ++j)
1645 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp2)(i, j) - (*tmp3)(i, j)) < tol, true);
1646 for(unsigned int j = 0 ; j < i; ++j)
1647 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(j, i)) < tol, true);
1648 }
1649 // ... Sparse
1650 add(*tmp2, * tmp4, *res);
1651 for(unsigned int i = 0; i < res->size(0); ++i)
1652 {
1653 for(unsigned int j = i ; j < res->size(1); ++j)
1654 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j) - (*tmp2)(i, j)) < tol, true);
1655 for(unsigned int j = 0 ; j < i; ++j)
1656 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j)) < tol, true);
1657 }
1658
1659 // ... Banded
1660 add(*tmp2, * tmp5, *res);
1661 for(signed i = 0; i < signed(Band->size1()); ++ i)
1662 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1663 if(j >= i)
1664 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp2)(i, j) - (*Band)(i, j)) < tol, true);
1665 else
1666 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*Band)(i, j)) < tol, true);
1667
1668 // ... Zero
1669 add(*tmp2, * tmp6, *res);
1670 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", *res == *tmp2, true);
1671
1672 // ... Identity
1673 add(*tmp2, * tmp7, *res);
1674 for(unsigned int i = 0; i < res->size(0); ++i)
1675 {
1676 for(unsigned int j = i ; j < res->size(1); ++j)
1677 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp2)(i, j) - (*tmp7)(i, j)) < tol, true);
1678 for(unsigned int j = 0 ; j < i; ++j)
1679 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j)) < tol, true);
1680 }
1681
1682 // triang - ...
1683 // ... dense
1684 sub(*tmp2, * tmp1, *res);
1685 for(unsigned int i = 0; i < res->size(0); ++i)
1686 {
1687 for(unsigned int j = i ; j < res->size(1); ++j)
1688 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp2)(i, j) + (*tmp1)(i, j)) < tol, true);
1689 for(unsigned int j = 0 ; j < i; ++j)
1690 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp1)(i, j)) < tol, true);
1691 }
1692 // ... Sym
1693 sub(*tmp2, * tmp3, *res);
1694 for(unsigned int i = 0; i < res->size(0); ++i)
1695 {
1696 for(unsigned int j = i ; j < res->size(1); ++j)
1697 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp2)(i, j) + (*tmp3)(i, j)) < tol, true);
1698 for(unsigned int j = 0 ; j < i; ++j)
1699 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp3)(j, i)) < tol, true);
1700 }
1701 // ... Sparse
1702 sub(*tmp2, * tmp4, *res);
1703 for(unsigned int i = 0; i < res->size(0); ++i)
1704 {
1705 for(unsigned int j = i ; j < res->size(1); ++j)
1706 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp2)(i, j) + (*tmp4)(i, j)) < tol, true);
1707 for(unsigned int j = 0 ; j < i; ++j)
1708 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp4)(i, j)) < tol, true);
1709 }
1710
1711 // ... Banded
1712 sub(*tmp2, * tmp5, *res);
1713 for(signed i = 0; i < signed(Band->size1()); ++ i)
1714 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1715 if(j >= i)
1716 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp2)(i, j) + (*Band)(i, j)) < tol, true);
1717 else
1718 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*Band)(i, j)) < tol, true);
1719
1720 // ... Zero
1721 sub(*tmp2, * tmp6, *res);
1722 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", *res == *tmp2, true);
1723
1724 // Identity
1725 sub(*tmp2, * tmp7, *res);
1726 for(unsigned int i = 0; i < res->size(0); ++i)
1727 {
1728 for(unsigned int j = i ; j < res->size(1); ++j)
1729 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp2)(i, j) + (*tmp7)(i, j)) < tol, true);
1730 for(unsigned int j = 0 ; j < i; ++j)
1731 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j)) < tol, true);
1732 }
1733
1734 // sym + ...
1735 // ... dense
1736 add(*tmp3, * tmp1, *res);
1737 for(unsigned int i = 0; i < res->size(0); ++i)
1738 {
1739 for(unsigned int j = i ; j < res->size(1); ++j)
1740 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) - (*tmp3)(i, j)) < tol, true);
1741 for(unsigned int j = 0 ; j < i; ++j)
1742 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) - (*tmp3)(j, i)) < tol, true);
1743 }
1744 // ... triang
1745 add(*tmp3, * tmp2, *res);
1746 for(unsigned int i = 0; i < res->size(0); ++i)
1747 {
1748 for(unsigned int j = i ; j < res->size(1); ++j)
1749 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(i, j) - (*tmp2)(i, j)) < tol, true);
1750 for(unsigned int j = 0 ; j < i; ++j)
1751 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(j, i)) < tol, true);
1752 }
1753 // ... Sparse
1754 add(*tmp3, * tmp4, *res);
1755 for(unsigned int i = 0; i < res->size(0); ++i)
1756 {
1757 for(unsigned int j = i ; j < res->size(1); ++j)
1758 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j) - (*tmp3)(i, j)) < tol, true);
1759 for(unsigned int j = 0 ; j < i; ++j)
1760 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j) - (*tmp3)(j, i)) < tol, true);
1761 }
1762
1763 // ... Banded
1764 add(*tmp3, * tmp5, *res);
1765 for(signed i = 0; i < signed(Band->size1()); ++ i)
1766 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1767 if(j >= i)
1768 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(i, j) - (*Band)(i, j)) < tol, true);
1769 else
1770 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(j, i) - (*Band)(i, j)) < tol, true);
1771
1772
1773 // ... Zero
1774 add(*tmp3, * tmp6, *res);
1775 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", *res == *tmp3, true);
1776
1777 // ... identity
1778 add(*tmp3, * tmp7, *res);
1779 for(unsigned int i = 0; i < res->size(0); ++i)
1780 {
1781 for(unsigned int j = i ; j < res->size(1); ++j)
1782 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp7)(i, j) - (*tmp3)(i, j)) < tol, true);
1783 for(unsigned int j = 0 ; j < i; ++j)
1784 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp7)(i, j) - (*tmp3)(j, i)) < tol, true);
1785 }
1786
1787 // sym - ...
1788 // ... dense
1789 sub(*tmp3, * tmp1, *res);
1790 for(unsigned int i = 0; i < res->size(0); ++i)
1791 {
1792 for(unsigned int j = i ; j < res->size(1); ++j)
1793 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(i, j) + (*tmp1)(i, j)) < tol, true);
1794 for(unsigned int j = 0 ; j < i; ++j)
1795 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(j, i) + (*tmp1)(i, j)) < tol, true);
1796 }
1797 // ... triang
1798 sub(*tmp3, * tmp2, *res);
1799 for(unsigned int i = 0; i < res->size(0); ++i)
1800 {
1801 for(unsigned int j = i ; j < res->size(1); ++j)
1802 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(i, j) + (*tmp2)(i, j)) < tol, true);
1803 for(unsigned int j = 0 ; j < i; ++j)
1804 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(j, i)) < tol, true);
1805 }
1806 // ... Sparse
1807 sub(*tmp3, * tmp4, *res);
1808 for(unsigned int i = 0; i < res->size(0); ++i)
1809 {
1810 for(unsigned int j = i ; j < res->size(1); ++j)
1811 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(i, j) + (*tmp4)(i, j)) < tol, true);
1812 for(unsigned int j = 0 ; j < i; ++j)
1813 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(j, i) + (*tmp4)(i, j)) < tol, true);
1814 }
1815
1816 // ... Banded
1817 sub(*tmp3, * tmp5, *res);
1818 for(signed i = 0; i < signed(Band->size1()); ++ i)
1819 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1820 if(j >= i)
1821 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(i, j) + (*Band)(i, j)) < tol, true);
1822 else
1823 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(j, i) + (*Band)(i, j)) < tol, true);
1824
1825 // ... Zero
1826 sub(*tmp3, * tmp6, *res);
1827 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", *res == *tmp3, true);
1828 // Identity
1829 sub(*tmp3, * tmp7, *res);
1830 for(unsigned int i = 0; i < res->size(0); ++i)
1831 {
1832 for(unsigned int j = i ; j < res->size(1); ++j)
1833 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(i, j) + (*tmp7)(i, j)) < tol, true);
1834 for(unsigned int j = 0 ; j < i; ++j)
1835 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(j, i) + (*tmp7)(i, j)) < tol, true);
1836 }
1837
1838 // sparse + ...
1839 // ... dense
1840 add(*tmp4, * tmp1, *res);
1841 for(unsigned int i = 0; i < res->size(0); ++i)
1842 for(unsigned int j = 0 ; j < res->size(1); ++j)
1843 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) - (*tmp4)(i, j)) < tol, true);
1844 // ... triang
1845 add(*tmp4, * tmp2, *res);
1846 for(unsigned int i = 0; i < res->size(0); ++i)
1847 {
1848 for(unsigned int j = i ; j < res->size(1); ++j)
1849 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j) - (*tmp2)(i, j)) < tol, true);
1850 for(unsigned int j = 0 ; j < i; ++j)
1851 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j)) < tol, true);
1852 }
1853 // ... Sym
1854 add(*tmp4, * tmp3, *res);
1855 for(unsigned int i = 0; i < res->size(0); ++i)
1856 {
1857 for(unsigned int j = i ; j < res->size(1); ++j)
1858 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j) - (*tmp3)(i, j)) < tol, true);
1859 for(unsigned int j = 0 ; j < i; ++j)
1860 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j) - (*tmp3)(j, i)) < tol, true);
1861 }
1862 // ... Banded
1863 add(*tmp4, * tmp5, *res);
1864 for(signed i = 0; i < signed(Band->size1()); ++ i)
1865 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1866 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j) - (*Band)(i, j)) < tol, true);
1867
1868 // ... zero
1869 add(*tmp4, * tmp6, *res);
1870 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", *res == *tmp4, true);
1871
1872 // sparse - ...
1873 // ... dense
1874 sub(*tmp4, * tmp1, *res);
1875 for(unsigned int i = 0; i < res->size(0); ++i)
1876 for(unsigned int j = 0 ; j < res->size(1); ++j)
1877 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j) + (*tmp1)(i, j)) < tol, true);
1878 // ... triangular
1879 sub(*tmp4, * tmp2, *res);
1880 for(unsigned int i = 0; i < res->size(0); ++i)
1881 {
1882 for(unsigned int j = i ; j < res->size(1); ++j)
1883 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j) + (*tmp2)(i, j)) < tol, true);
1884 for(unsigned int j = 0 ; j < i; ++j)
1885 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j)) < tol, true);
1886 }
1887 // ... Sym
1888 sub(*tmp4, * tmp3, *res);
1889 for(unsigned int i = 0; i < res->size(0); ++i)
1890 {
1891 for(unsigned int j = i ; j < res->size(1); ++j)
1892 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j) + (*tmp3)(i, j)) < tol, true);
1893 for(unsigned int j = 0 ; j < i; ++j)
1894 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j) + (*tmp3)(j, i)) < tol, true);
1895 }
1896
1897 // ... Banded
1898 sub(*tmp4, * tmp5, *res);
1899 for(signed i = 0; i < signed(Band->size1()); ++ i)
1900 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1901 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j) + (*Band)(i, j)) < tol, true);
1902
1903 // ... zero
1904 sub(*tmp4, * tmp6, *res);
1905 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", *res == *tmp4, true);
1906
1907 // Banded + ...
1908 // ... dense
1909 add(*tmp5, * tmp1, *res);
1910 for(signed i = 0; i < signed(Band->size1()); ++ i)
1911 {
1912 for(signed j = 0; j < std::max(i - 1, 0); ++j)
1913 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j)) < tol, true);
1914 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1915 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j) - (*tmp5)(i, j)) < tol, true);
1916 for(signed j = std::min(i + 2, signed(Band->size2())); j < signed(Band->size2()); ++j)
1917 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp1)(i, j)) < tol, true);
1918 }
1919 // ... triang
1920 add(*tmp5, * tmp2, *res);
1921 for(signed i = 0; i < signed(Band->size1()); ++ i)
1922 {
1923 for(signed j = 0; j < std::max(i - 1, 0); ++j)
1924 if(j >= i) CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp2)(i, j)) < tol, true);
1925 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1926 if(j >= i)
1927 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp2)(i, j) - (*tmp5)(i, j)) < tol, true);
1928 else
1929 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp5)(i, j)) < tol, true);
1930 for(signed j = std::min(i + 2, signed(Band->size2())); j < signed(Band->size2()); ++j)
1931 if(j >= i) CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp2)(i, j)) < tol, true);
1932 }
1933
1934 // ...sym
1935 add(*tmp5, * tmp3, *res);
1936 for(signed i = 0; i < signed(Band->size1()); ++ i)
1937 {
1938 for(signed j = 0; j < std::max(i - 1, 0); ++j)
1939 if(j >= i) CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(i, j)) < tol, true);
1940 else CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(j, i)) < tol, true);
1941 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1942 if(j >= i)
1943 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(i, j) - (*tmp5)(i, j)) < tol, true);
1944 else
1945 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(j, i) - (*tmp5)(i, j)) < tol, true);
1946 for(signed j = std::min(i + 2, signed(Band->size2())); j < signed(Band->size2()); ++j)
1947 if(j >= i) CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(i, j)) < tol, true);
1948 else CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp3)(j, i)) < tol, true);
1949 }
1950
1951 //... sparse
1952 add(*tmp5, * tmp4, *res);
1953 for(signed i = 0; i < signed(Band->size1()); ++ i)
1954 {
1955 for(signed j = 0; j < std::max(i - 1, 0); ++j)
1956 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j)) < tol, true);
1957 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1958 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j) - (*tmp5)(i, j)) < tol, true);
1959 for(signed j = std::min(i + 2, signed(Band->size1())); j < signed(Band->size1()); ++j)
1960 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp4)(i, j)) < tol, true);
1961 }
1962
1963 // ... zero
1964 add(*tmp5, * tmp6, *res);
1965 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", *res == *tmp5, true);
1966 // ... identity
1967 add(*tmp5, * tmp7, *res);
1968 for(signed i = 0; i < signed(Band->size1()); ++ i)
1969 {
1970 for(signed j = 0; j < std::max(i - 1, 0); ++j)
1971 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp7)(i, j)) < tol, true);
1972 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1973 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp7)(i, j) - (*tmp5)(i, j)) < tol, true);
1974 for(signed j = std::min(i + 2, signed(Band->size2())); j < signed(Band->size2()); ++j)
1975 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp7)(i, j)) < tol, true);
1976 }
1977
1978 // Banded - ...
1979 // ... dense
1980
1981 sub(*tmp5, * tmp1, *res);
1982 for(signed i = 0; i < signed(Band->size1()); ++ i)
1983 {
1984 for(signed j = 0; j < std::max(i - 1, 0); ++j)
1985 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp1)(i, j)) < tol, true);
1986 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1987 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp5)(i, j) + (*tmp1)(i, j)) < tol, true);
1988 for(signed j = std::min(i + 2, signed(Band->size2())); j < signed(Band->size2()); ++j)
1989 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp1)(i, j)) < tol, true);
1990 }
1991
1992 // ... triang
1993 sub(*tmp5, * tmp2, *res);
1994 for(signed i = 0; i < signed(Band->size1()); ++ i)
1995 {
1996 for(signed j = 0; j < std::max(i - 1, 0); ++j)
1997 if(j >= i) CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp2)(i, j)) < tol, true);
1998 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
1999 if(j >= i)
2000 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp5)(i, j) + (*tmp2)(i, j)) < tol, true);
2001 else
2002 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp5)(i, j)) < tol, true);
2003 for(signed j = std::min(i + 2, signed(Band->size2())); j < signed(Band->size2()); ++j)
2004 if(j >= i) CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp2)(i, j)) < tol, true);
2005 }
2006
2007 // ...sym
2008 sub(*tmp5, * tmp3, *res);
2009 for(signed i = 0; i < signed(Band->size1()); ++ i)
2010 {
2011 for(signed j = 0; j < std::max(i - 1, 0); ++j)
2012 if(j >= i) CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp3)(i, j)) < tol, true);
2013 else CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp3)(j, i)) < tol, true);
2014 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
2015 if(j >= i)
2016 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp5)(i, j) + (*tmp3)(i, j)) < tol, true);
2017 else
2018 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp5)(i, j) + (*tmp3)(j, i)) < tol, true);
2019 for(signed j = std::min(i + 2, signed(Band->size2())); j < signed(Band->size2()); ++j)
2020 if(j >= i) CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp3)(i, j)) < tol, true);
2021 else CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp3)(j, i)) < tol, true);
2022 }
2023
2024 //... sparse
2025 sub(*tmp5, * tmp4, *res);
2026 for(signed i = 0; i < signed(Band->size1()); ++ i)
2027 {
2028 for(signed j = 0; j < std::max(i - 1, 0); ++j)
2029 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp4)(i, j)) < tol, true);
2030 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
2031 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp4)(i, j) - (*tmp5)(i, j)) < tol, true);
2032 for(signed j = std::min(i + 2, signed(Band->size1())); j < signed(Band->size1()); ++j)
2033 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp4)(i, j)) < tol, true);
2034 }
2035
2036 // ... zero
2037 sub(*tmp5, * tmp6, *res);
2038 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", *res == *tmp5, true);
2039 // ... identity
2040 sub(*tmp5, * tmp7, *res);
2041 for(signed i = 0; i < signed(Band->size1()); ++ i)
2042 {
2043 for(signed j = 0; j < std::max(i - 1, 0); ++j)
2044 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp7)(i, j)) < tol, true);
2045 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(Band->size2())); ++ j)
2046 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) - (*tmp5)(i, j) + (*tmp7)(i, j)) < tol, true);
2047 for(signed j = std::min(i + 2, signed(Band->size2())); j < signed(Band->size2()); ++j)
2048 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*res)(i, j) + (*tmp7)(i, j)) < tol, true);
2049 }
2050
2051 std::cout << "--> test operators7 ended with success." <<std::endl;
2052 }
2053
2054
2055
2056 //void SimpleMatrixTest::testOperators8()
2057 // {
2058 // std::cout << "--> Test: operator8." <<std::endl;
2059
2060 // // // Simple = Simple * Simple
2061 // // *C = prod(*A, *B);
2062 // // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8: ", norm_inf(*C->dense() - prod(*A->dense(), *B->dense())) < tol, true);
2063
2064 // // Block = Simple * Simple
2065 // // *Cb = prod(*A, *B);
2066 // // DenseMat Dtmp = prod(*A->dense(), *B->dense());
2067 // // SP::SimpleMatrix tmp(new SimpleMatrix(Dtmp));
2068 // // for (unsigned int i = 0; i < C->size(0); ++i)
2069 // // for (unsigned int j = i ; j < C->size(1); ++j)
2070 // // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*Cb)(i, j) - (*tmp)(i, j)) < tol, true);
2071
2072 // // Others ...
2073
2074 // SP::SiconosMatrix tmp1(new SimpleMatrix(4, 4, 2.3));
2075 // SP::SiconosMatrix tmp2(new SimpleMatrix(*T2));
2076 // SP::SiconosMatrix tmp3(new SimpleMatrix(*S2));
2077 // SP::SiconosMatrix tmp4(new SimpleMatrix(*SP));
2078 // SP::SiconosMatrix tmp5(new SimpleMatrix(*Band));
2079
2080 // SP::SiconosMatrix res(new SimpleMatrix(4, 4, 0));
2081
2082 // // Dense * ...
2083 // // triang
2084 // *res = prod(*tmp1, *tmp2);
2085 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp1->getDense(), tmp2->getTriang())) < tol, true);
2086 // // Sym
2087 // *res = prod(*tmp1, *tmp3);
2088 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp1->getDense(), tmp3->getSym())) < tol, true);
2089 // // Sparse
2090 // *res = prod(*tmp1, *tmp4);
2091 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp1->getDense(), tmp4->getSparse())) < tol, true);
2092 // // Banded
2093 // *res = prod(*tmp1, *tmp5);
2094 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp1->getDense(), tmp5->getBanded())) < tol, true);
2095 // // triang * ...
2096 // // dense
2097 // *res = prod(*tmp2, *tmp1);
2098 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp2->getTriang(), tmp1->getDense())) < tol, true);
2099 // // Sym
2100 // *res = prod(*tmp2, *tmp3);
2101 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp2->getTriang(), tmp3->getSym())) < tol, true);
2102 // // Sparse
2103 // *res = prod(*tmp2, *tmp4);
2104 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp2->getTriang(), tmp4->getSparse())) < tol, true);
2105 // // Banded
2106 // *res = prod(*tmp2, *tmp5);
2107 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp2->getTriang(), tmp5->getBanded())) < tol, true);
2108 // // sym * ...
2109 // // dense
2110 // *res = prod(*tmp3, *tmp1);
2111 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp3->getSym(), tmp1->getDense())) < tol, true);
2112 // // triang
2113 // *res = prod(*tmp3, *tmp2);
2114 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp3->getSym(), tmp2->getTriang())) < tol, true);
2115 // // Sparse
2116 // *res = prod(*tmp3, *tmp4);
2117 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp3->getSym(), tmp4->getSparse())) < tol, true);
2118 // // Banded
2119 // *res = prod(*tmp3, *tmp5);
2120 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp3->getSym(), tmp5->getBanded())) < tol, true);
2121 // // Sparse * ...
2122 // // dense
2123 // *res = prod(*tmp4, *tmp1);
2124 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp4->getSparse(), tmp1->getDense())) < tol, true);
2125 // // triang
2126 // *res = prod(*tmp4, *tmp2);
2127 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp4->getSparse(), tmp2->getTriang())) < tol, true);
2128 // // Sym
2129 // *res = prod(*tmp4, *tmp3);
2130 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp4->getSparse(), tmp3->getSym())) < tol, true);
2131 // // Banded
2132 // *res = prod(*tmp4, *tmp5);
2133 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp4->getSparse(), tmp5->getBanded())) < tol, true);
2134 // // Banded * ...
2135 // // dense
2136 // *res = prod(*tmp5, *tmp1);
2137 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp5->getBanded(), tmp1->getDense())) < tol, true);
2138 // // triang
2139 // *res = prod(*tmp5, *tmp2);
2140 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp5->getBanded(), tmp2->getTriang())) < tol, true);
2141 // // Sparse
2142 // *res = prod(*tmp5, *tmp4);
2143 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp5->getBanded(), tmp4->getSparse())) < tol, true);
2144 // // Sym
2145 // *res = prod(*tmp5, *tmp3);
2146 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(*res->dense() - prod(tmp5->getBanded(), tmp3->getSym())) < tol, true);
2147
2148 // std::cout << "--> test operators8 ended with success." <<std::endl;
2149 // }
2150
testOperators8Bis()2151 void SimpleMatrixTest::testOperators8Bis()
2152 {
2153 std::cout << "--> Test: operator8Bis." <<std::endl;
2154 // Simple = Simple * Simple
2155 prod(*A, *B, *C);
2156 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*C->dense() - prod(*A->dense(), *B->dense())) < tol, true);
2157
2158 // Block = Simple * Simple
2159 prod(*A, *B, *Cb);
2160 DenseMat Dtmp = prod(*A->dense(), *B->dense());
2161 SP::SimpleMatrix tmp(new SimpleMatrix(Dtmp));
2162 for(unsigned int i = 0; i < Cb->size(0); ++i)
2163 for(unsigned int j = i ; j < Cb->size(1); ++j)
2164 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*Cb)(i, j) - (*tmp)(i, j)) < tol, true);
2165
2166 // Others ...
2167
2168 // Others ...
2169 SP::SiconosMatrix tmp1(new SimpleMatrix(4, 4, 2.4));
2170 SP::SiconosMatrix tmp2(new SimpleMatrix(*T2));
2171 SP::SiconosMatrix tmp3(new SimpleMatrix(*S2));
2172 SP::SiconosMatrix tmp4(new SimpleMatrix(*SP));
2173 SP::SiconosMatrix tmp5(new SimpleMatrix(*Band));
2174
2175 SP::SiconosMatrix res(new SimpleMatrix(4, 4));
2176
2177 // Dense * ...
2178 // triang
2179 prod(*tmp1, *tmp2, *res);
2180 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp1->getDense(), tmp2->getTriang())) < tol, true);
2181 // Sym
2182 prod(*tmp1, *tmp3, *res);
2183 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp1->getDense(), tmp3->getSym())) < tol, true);
2184 // Sparse
2185 prod(*tmp1, *tmp4, *res);
2186 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp1->getDense(), tmp4->getSparse())) < tol, true);
2187 // Banded
2188 prod(*tmp1, *tmp5, *res);
2189 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp1->getDense(), tmp5->getBanded())) < tol, true);
2190 // triang * ...
2191 // dense
2192 prod(*tmp2, *tmp1, *res);
2193 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp2->getTriang(), tmp1->getDense())) < tol, true);
2194 // Sym
2195 prod(*tmp2, *tmp3, *res);
2196 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp2->getTriang(), tmp3->getSym())) < tol, true);
2197 // Sparse
2198 prod(*tmp2, *tmp4, *res);
2199 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp2->getTriang(), tmp4->getSparse())) < tol, true);
2200 // Banded
2201 prod(*tmp2, *tmp5, *res);
2202 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp2->getTriang(), tmp5->getBanded())) < tol, true);
2203 // sym * ...
2204 // dense
2205 prod(*tmp3, *tmp1, *res);
2206 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp3->getSym(), tmp1->getDense())) < tol, true);
2207 // triang
2208 prod(*tmp3, *tmp2, *res);
2209 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp3->getSym(), tmp2->getTriang())) < tol, true);
2210 // Sparse
2211 prod(*tmp3, *tmp4, *res);
2212 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp3->getSym(), tmp4->getSparse())) < tol, true);
2213 // Banded
2214 prod(*tmp3, *tmp5, *res);
2215 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp3->getSym(), tmp5->getBanded())) < tol, true);
2216 // Sparse * ...
2217 // dense
2218 prod(*tmp4, *tmp1, *res);
2219 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp4->getSparse(), tmp1->getDense())) < tol, true);
2220 // triang
2221 prod(*tmp4, *tmp2, *res);
2222 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp4->getSparse(), tmp2->getTriang())) < tol, true);
2223 // Sym
2224 prod(*tmp4, *tmp3, *res);
2225 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp4->getSparse(), tmp3->getSym())) < tol, true);
2226 // Banded
2227 prod(*tmp4, *tmp5, *res);
2228 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp4->getSparse(), tmp5->getBanded())) < tol, true);
2229 // Banded * ...
2230 // dense
2231 prod(*tmp5, *tmp1, *res);
2232 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp5->getBanded(), tmp1->getDense())) < tol, true);
2233 // triang
2234 prod(*tmp5, *tmp2, *res);
2235 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp5->getBanded(), tmp2->getTriang())) < tol, true);
2236 // Sparse
2237 prod(*tmp5, *tmp4, *res);
2238 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp5->getBanded(), tmp4->getSparse())) < tol, true);
2239 // Sym
2240 prod(*tmp5, *tmp3, *res);
2241 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Bis: ", norm_inf(*res->dense() - prod(tmp5->getBanded(), tmp3->getSym())) < tol, true);
2242
2243 std::cout << "--> test operators8Bis ended with success." <<std::endl;
2244 }
2245
testOperators8Ter()2246 void SimpleMatrixTest::testOperators8Ter()
2247 {
2248 std::cout << "--> Test: operator8Ter." <<std::endl;
2249 // Simple = Simple * Simple
2250 axpy_prod(*A, *B, *C, true);
2251 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Ter: ", norm_inf(*C->dense() - prod(*A->dense(), *B->dense())) < tol, true);
2252
2253 // Simple += Simple * Simple
2254 SP::SiconosMatrix backUp(new SimpleMatrix(*C));
2255
2256 axpy_prod(*A, *B, *C, false);
2257
2258 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8Ter: ", norm_inf(*C->dense() - prod(*A->dense(), *B->dense()) - *backUp->dense()) < tol, true);
2259 // Block = Simple * Simple
2260 axpy_prod(*A, *B, *Cb, true);
2261 DenseMat Dtmp = prod(*A->dense(), *B->dense());
2262 SP::SimpleMatrix tmp(new SimpleMatrix(Dtmp));
2263 for(unsigned int i = 0; i < Cb->size(0); ++i)
2264 for(unsigned int j = i ; j < Cb->size(1); ++j)
2265 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*Cb)(i, j) - (*tmp)(i, j)) < tol, true);
2266
2267 *backUp = *Cb;
2268 // Block += Simple * Simple
2269 axpy_prod(*A, *B, *Cb, false);
2270 Dtmp = prod(*A->dense(), *B->dense());
2271 *tmp = Dtmp;
2272 for(unsigned int i = 0; i < Cb->size(0); ++i)
2273 for(unsigned int j = i ; j < Cb->size(1); ++j)
2274 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*Cb)(i, j) - (*tmp)(i, j) - (*backUp)(i, j)) < tol, true);
2275
2276 std::cout << "--> test operators8Ter ended with success." <<std::endl;
2277 }
2278
testOperators8_4()2279 void SimpleMatrixTest::testOperators8_4() // C += A*B
2280 {
2281 std::cout << "--> Test: operator8_4." <<std::endl;
2282 // Simple = Simple * Simple
2283 C->zero();
2284 prod(*A, *B, *C, false);
2285 prod(*A, *B, *C, false);
2286 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_4: ", norm_inf(*C->dense() - 2 * prod(*A->dense(), *B->dense())) < tol, true);
2287
2288 // Block = Simple * Simple
2289 Cb->zero();
2290 prod(*A, *B, *Cb, false);
2291 prod(*A, *B, *Cb, false);
2292 DenseMat Dtmp = prod(*A->dense(), *B->dense());
2293 SP::SimpleMatrix tmp(new SimpleMatrix(Dtmp));
2294 for(unsigned int i = 0; i < Cb->size(0); ++i)
2295 for(unsigned int j = i ; j < Cb->size(1); ++j)
2296 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", fabs((*Cb)(i, j) - 2 * (*tmp)(i, j)) < tol, true);
2297 std::cout << "--> test operators8_4 ended with success." <<std::endl;
2298 }
2299
testOperators8_5()2300 void SimpleMatrixTest::testOperators8_5()
2301 {
2302 // == Test subprod ==
2303
2304 std::cout << "--> Test: operator8_5." <<std::endl;
2305 Index coord(8);
2306 SP::SiconosVector x1(new SiconosVector(2));
2307 SP::SiconosVector x2(new SiconosVector(3));
2308 SP::SiconosVector x3(new SiconosVector(5));
2309 SP::SiconosVector y(new SiconosVector(size));
2310 SP::BlockVector x(new BlockVector());
2311 SP::SiconosVector v(new SiconosVector(size));
2312 x->insertPtr(x1);
2313 x->insertPtr(x2);
2314 x->insertPtr(x3);
2315 for(unsigned int i = 0 ; i < size; ++i)
2316 {
2317 (*x)(i) = (double)i + 3;
2318 (*v)(i) = (double)i + 3;
2319 }
2320
2321 // v == x but x is a 3-blocks vector.
2322
2323 // Simple = Simple * Simple, all dense
2324 // subprod but with full matrix/vectors
2325 coord[0] = 0;
2326 coord[1] = size;
2327 coord[2] = 0;
2328 coord[3] = size;
2329 coord[4] = 0;
2330 coord[5] = size;
2331 coord[6] = 0;
2332 coord[7] = size;
2333 subprod(*A, *v, *y, coord, true);
2334 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", norm_inf(*y->dense() - prod(*A->dense(), *v->dense())) < tol, true);
2335
2336 // Simple = Simple * Block, all dense
2337 // subprod but with full matrix/vectors
2338 // subprod(*A,*x,*y, coord, true);
2339 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", norm_inf(*y->dense()- prod(*A->dense(),*v->dense()))<tol, true);
2340
2341 coord[0] = 0;
2342 coord[1] = 2;
2343 coord[2] = 1;
2344 coord[3] = 3;
2345 coord[4] = 3;
2346 coord[5] = 5;
2347 coord[6] = 2;
2348 coord[7] = 4;
2349 y->zero();
2350 // Simple = Simple * Simple, all dense
2351 subprod(*A, *v, *y, coord, true);
2352 double res = (*A)(0, 1) * (*v)(3) + (*A)(0, 2) * (*v)(4);
2353 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs(res - (*y)(2)) < tol, true);
2354 res = (*A)(1, 1) * (*v)(3) + (*A)(1, 2) * (*v)(4);
2355 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs(res - (*y)(3)) < tol, true);
2356 for(unsigned int i = 0; i < size; ++i)
2357 {
2358 if(i != 2 && i != 3)
2359 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs((*y)(i)) < tol, true);
2360 }
2361 y->zero();
2362 // Simple = Simple * Block, all dense
2363 // subprod(*A,*x,*y, coord, true);
2364 // res = (*A)(0,1)*(*x)(3) + (*A)(0,2)*(*x)(4);
2365 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs(res-(*y)(2))<tol, true);
2366 // res = (*A)(1,1)*(*x)(3) + (*A)(1,2)*(*x)(4);
2367 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs(res-(*y)(3))<tol, true);
2368 // for (unsigned int i=0; i<size; ++i)
2369 // {
2370 // if (i!=2 && i!=3)
2371 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs((*y)(i))<tol, true);
2372 // }
2373 // // Others ...
2374 // Triang
2375
2376 SP::SiconosMatrix A2(new SimpleMatrix(10, 10, TRIANGULAR));
2377 for(unsigned i = 0; i < A2->size(0); ++ i)
2378 for(unsigned j = i; j < A2->size(1); ++ j)
2379 (*A2)(i, j) = 3 * i + j;
2380
2381 subprod(*A2, *v, *y, coord, true);
2382 res = (*A2)(0, 1) * (*v)(3) + (*A2)(0, 2) * (*v)(4);
2383 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs(res - (*y)(2)) < tol, true);
2384 res = (*A2)(1, 1) * (*v)(3) + (*A2)(1, 2) * (*v)(4);
2385 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs(res - (*y)(3)) < tol, true);
2386 for(unsigned int i = 0; i < size; ++i)
2387 {
2388 if(i != 2 && i != 3)
2389 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs((*y)(i)) < tol, true);
2390 }
2391 // Sym
2392 A2.reset(new SimpleMatrix(10, 10, SYMMETRIC));
2393 for(unsigned i = 0; i < A2->size(0); ++ i)
2394 for(unsigned j = i; j < A2->size(1); ++ j)
2395 (*A2)(i, j) = 3 * i + j;
2396
2397 subprod(*A2, *v, *y, coord, true);
2398 res = (*A2)(0, 1) * (*v)(3) + (*A2)(0, 2) * (*v)(4);
2399 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs(res - (*y)(2)) < tol, true);
2400 res = (*A2)(1, 1) * (*v)(3) + (*A2)(1, 2) * (*v)(4);
2401 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs(res - (*y)(3)) < tol, true);
2402 for(unsigned int i = 0; i < size; ++i)
2403 {
2404 if(i != 2 && i != 3)
2405 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs((*y)(i)) < tol, true);
2406 }
2407
2408 // Sparse
2409 A2.reset(new SimpleMatrix(10, 10, Siconos::SPARSE));
2410 for(unsigned i = 0; i < A2->size(0); ++ i)
2411 for(unsigned j = i; j < A2->size(1); ++ j)
2412 A2->setValue(i, j, 3 * i + j);
2413
2414 subprod(*A2, *v, *y, coord, true);
2415 res = (*A2)(0, 1) * (*v)(3) + (*A2)(0, 2) * (*v)(4);
2416 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs(res - (*y)(2)) < tol, true);
2417 res = (*A2)(1, 1) * (*v)(3) + (*A2)(1, 2) * (*v)(4);
2418 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs(res - (*y)(3)) < tol, true);
2419 for(unsigned int i = 0; i < size; ++i)
2420 {
2421 if(i != 2 && i != 3)
2422 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs((*y)(i)) < tol, true);
2423 }
2424
2425 // Banded
2426 A2.reset(new SimpleMatrix(10, 10, BANDED));
2427 for(signed i = 0; i < signed(A2->size(0)); ++ i)
2428 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(A2->size(1))); ++ j)
2429 (*A2)(i, j) = 3 * i + j;
2430 subprod(*A2, *v, *y, coord, true);
2431 res = (*A2)(0, 1) * (*v)(3);
2432 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs(res - (*y)(2)) < tol, true);
2433 res = (*A2)(1, 1) * (*v)(3) + (*A2)(1, 2) * (*v)(4);
2434 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs(res - (*y)(3)) < tol, true);
2435 for(unsigned int i = 0; i < size; ++i)
2436 {
2437 if(i != 2 && i != 3)
2438 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_5: ", fabs((*y)(i)) < tol, true);
2439 }
2440
2441 std::cout << "--> test operators8_5 ended with success." <<std::endl;
2442 }
2443
testOperators8_6()2444 void SimpleMatrixTest::testOperators8_6()
2445 {
2446 // == Test subprod, with += ==
2447
2448 std::cout << "--> Test: operator8_6." <<std::endl;
2449 Index coord(8);
2450 SP::SiconosVector x1(new SiconosVector(2));
2451 SP::SiconosVector x2(new SiconosVector(3));
2452 SP::SiconosVector x3(new SiconosVector(5));
2453 SP::SiconosVector y(new SiconosVector(size));
2454 SP::BlockVector x(new BlockVector());
2455 SP::SiconosVector v(new SiconosVector(size));
2456 x->insertPtr(x1);
2457 x->insertPtr(x2);
2458 x->insertPtr(x3);
2459 for(unsigned int i = 0 ; i < size; ++i)
2460 {
2461 (*x)(i) = (double)i + 3;
2462 (*v)(i) = (double)i + 3;
2463 }
2464
2465 // v == x but x is a 3-blocks vector.
2466
2467 *y = *v;
2468
2469 // Simple = Simple * Simple, all dense
2470 // subprod but with full matrix/vectors
2471 coord[0] = 0;
2472 coord[1] = size;
2473 coord[2] = 0;
2474 coord[3] = size;
2475 coord[4] = 0;
2476 coord[5] = size;
2477 coord[6] = 0;
2478 coord[7] = size;
2479 subprod(*A, *v, *y, coord, false);
2480 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", norm_inf(*y->dense() - prod(*A->dense(), *v->dense()) - *v->dense()) < tol, true);
2481
2482 // Simple = Simple * Block, all dense
2483 // subprod but with full matrix/vectors
2484 *y = *v;
2485 // subprod(*A,*x,*y, coord, false);
2486 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", norm_inf(*y->dense()- prod(*A->dense(),*v->dense())- *v->dense())<tol, true);
2487
2488 coord[0] = 0;
2489 coord[1] = 2;
2490 coord[2] = 1;
2491 coord[3] = 3;
2492 coord[4] = 3;
2493 coord[5] = 5;
2494 coord[6] = 2;
2495 coord[7] = 4;
2496
2497 // Simple = Simple * Simple, all dense
2498 *y = *v;
2499 subprod(*A, *v, *y, coord, false);
2500 double res = (*A)(0, 1) * (*v)(3) + (*A)(0, 2) * (*v)(4) + (*v)(2);
2501 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs(res - (*y)(2)) < tol, true);
2502 res = (*A)(1, 1) * (*v)(3) + (*A)(1, 2) * (*v)(4) + (*v)(3);
2503 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs(res - (*y)(3)) < tol, true);
2504 for(unsigned int i = 0; i < size; ++i)
2505 {
2506 if(i != 2 && i != 3)
2507 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs((*y)(i) - (*v)(i)) < tol, true);
2508 }
2509 *y = *v;
2510 // Simple = Simple * Block, all dense
2511 // subprod(*A,*x,*y, coord, false);
2512 // res = (*A)(0,1)*(*x)(3) + (*A)(0,2)*(*x)(4) + (*v)(2);
2513 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs(res-(*y)(2))<tol, true);
2514 // res = (*A)(1,1)*(*x)(3) + (*A)(1,2)*(*x)(4) + (*v)(3);
2515 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs(res-(*y)(3))<tol, true);
2516 // for (unsigned int i=0; i<size; ++i)
2517 // {
2518 // if (i!=2 && i!=3)
2519 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs((*y)(i)-(*v)(i))<tol, true);
2520 // }
2521
2522 // // Others ...
2523 // Triang
2524
2525 SP::SiconosMatrix A2(new SimpleMatrix(10, 10, TRIANGULAR));
2526 for(unsigned i = 0; i < A2->size(0); ++ i)
2527 for(unsigned j = i; j < A2->size(1); ++ j)
2528 (*A2)(i, j) = 3 * i + j;
2529
2530 *y = *v;
2531 subprod(*A2, *v, *y, coord, false);
2532 res = (*A2)(0, 1) * (*v)(3) + (*A2)(0, 2) * (*v)(4) + (*v)(2);
2533 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs(res - (*y)(2)) < tol, true);
2534 res = (*A2)(1, 1) * (*v)(3) + (*A2)(1, 2) * (*v)(4) + (*v)(3);
2535 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs(res - (*y)(3)) < tol, true);
2536 for(unsigned int i = 0; i < size; ++i)
2537 {
2538 if(i != 2 && i != 3)
2539 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs((*y)(i) - (*v)(i)) < tol, true);
2540 }
2541
2542 // Sym
2543 A2.reset(new SimpleMatrix(10, 10, SYMMETRIC));
2544 for(unsigned i = 0; i < A2->size(0); ++ i)
2545 for(unsigned j = i; j < A2->size(1); ++ j)
2546 (*A2)(i, j) = 3 * i + j;
2547
2548 *y = *v;
2549 subprod(*A2, *v, *y, coord, false);
2550 res = (*A2)(0, 1) * (*v)(3) + (*A2)(0, 2) * (*v)(4) + (*v)(2);
2551 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs(res - (*y)(2)) < tol, true);
2552 res = (*A2)(1, 1) * (*v)(3) + (*A2)(1, 2) * (*v)(4) + (*v)(3);
2553 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs(res - (*y)(3)) < tol, true);
2554 for(unsigned int i = 0; i < size; ++i)
2555 {
2556 if(i != 2 && i != 3)
2557 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs((*y)(i) - (*v)(i)) < tol, true);
2558 }
2559
2560 // Sparse
2561 A2.reset(new SimpleMatrix(10, 10, Siconos::SPARSE));
2562 for(unsigned i = 0; i < A2->size(0); ++ i)
2563 for(unsigned j = i; j < A2->size(1); ++ j)
2564 A2->setValue(i, j, 3 * i + j);
2565
2566 *y = *v;
2567 subprod(*A2, *v, *y, coord, false);
2568 res = (*A2)(0, 1) * (*v)(3) + (*A2)(0, 2) * (*v)(4) + (*v)(2);
2569 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs(res - (*y)(2)) < tol, true);
2570 res = (*A2)(1, 1) * (*v)(3) + (*A2)(1, 2) * (*v)(4) + (*v)(3);
2571 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs(res - (*y)(3)) < tol, true);
2572 for(unsigned int i = 0; i < size; ++i)
2573 {
2574 if(i != 2 && i != 3)
2575 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs((*y)(i) - (*v)(i)) < tol, true);
2576 }
2577
2578 // Banded
2579 A2.reset(new SimpleMatrix(10, 10, BANDED));
2580 for(signed i = 0; i < signed(A2->size(0)); ++ i)
2581 for(signed j = std::max(i - 1, 0); j < std::min(i + 2, signed(A2->size(1))); ++ j)
2582 (*A2)(i, j) = 3 * i + j;
2583 *y = *v;
2584 subprod(*A2, *v, *y, coord, false);
2585 res = (*A2)(0, 1) * (*v)(3) + (*v)(2);
2586 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs(res - (*y)(2)) < tol, true);
2587 res = (*A2)(1, 1) * (*v)(3) + (*A2)(1, 2) * (*v)(4) + (*v)(3);
2588 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs(res - (*y)(3)) < tol, true);
2589 for(unsigned int i = 0; i < size; ++i)
2590 {
2591 if(i != 2 && i != 3)
2592 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators8_6: ", fabs((*y)(i) - (*v)(i)) < tol, true);
2593 }
2594
2595 std::cout << "--> test operators8_6 ended with success." <<std::endl;
2596 }
2597
testOperators9()2598 void SimpleMatrixTest::testOperators9()
2599 {
2600 std::cout << "--> Test: operator9." <<std::endl;
2601
2602 // C = a*A or A/a
2603
2604 double a = 2.2;
2605 int a1 = 3;
2606
2607 // Simple = a * Simple or Simple/a
2608 *C = a * *A;
2609 for(unsigned int i = 0; i < C->size(0); ++i)
2610 for(unsigned int j = i ; j < C->size(1); ++j)
2611 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*C)(i, j) - a * (*A)(i, j)) < tol, true);
2612 *C = a1 * *A;
2613 for(unsigned int i = 0; i < C->size(0); ++i)
2614 for(unsigned int j = i ; j < C->size(1); ++j)
2615 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*C)(i, j) - a1 * (*A)(i, j)) < tol, true);
2616
2617 // *C = *A / a;
2618 // for (unsigned int i = 0; i < C->size(0); ++i)
2619 // for (unsigned int j = i ; j < C->size(1); ++j)
2620 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*C)(i, j) - (*A)(i, j) / a) < tol, true);
2621 // *C = *A / a1;
2622 // for (unsigned int i = 0; i < C->size(0); ++i)
2623 // for (unsigned int j = i ; j < C->size(1); ++j)
2624 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*C)(i, j) - (*A)(i, j) / a1) < tol, true);
2625
2626 // Simple = a * Block
2627
2628 *C = a * *Ab;
2629 for(unsigned int i = 0; i < C->size(0); ++i)
2630 for(unsigned int j = i ; j < C->size(1); ++j)
2631 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*C)(i, j) - a * (*Ab)(i, j)) < tol, true);
2632 ;
2633 *C = a1 * *Ab;
2634 for(unsigned int i = 0; i < C->size(0); ++i)
2635 for(unsigned int j = i ; j < C->size(1); ++j)
2636 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*C)(i, j) - a1 * (*Ab)(i, j)) < tol, true);
2637
2638 // *C = *Ab / a;
2639 // for (unsigned int i = 0; i < C->size(0); ++i)
2640 // for (unsigned int j = i ; j < C->size(1); ++j)
2641 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*C)(i, j) - (*Ab)(i, j) / a) < tol, true);
2642 // *C = *Ab / a1;
2643 // for (unsigned int i = 0; i < C->size(0); ++i)
2644 // for (unsigned int j = i ; j < C->size(1); ++j)
2645 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*C)(i, j) - (*Ab)(i, j) / a1) < tol, true);
2646
2647 // Block = a * Block
2648 *Cb = a * *Ab;
2649 for(unsigned int i = 0; i < Cb->size(0); ++i)
2650 for(unsigned int j = i ; j < Cb->size(1); ++j)
2651 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*Cb)(i, j) - a * (*Ab)(i, j)) < tol, true);
2652 *Cb = a1 * *Ab;
2653 for(unsigned int i = 0; i < Cb->size(0); ++i)
2654 for(unsigned int j = i ; j < Cb->size(1); ++j)
2655 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*Cb)(i, j) - a1 * (*Ab)(i, j)) < tol, true);
2656
2657 // *Cb = *Ab / a;
2658 // for (unsigned int i = 0; i < Cb->size(0); ++i)
2659 // for (unsigned int j = i ; j < Cb->size(1); ++j)
2660 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*Cb)(i, j) - (*Ab)(i, j) / a) < tol, true);
2661 // *Cb = *Ab / a1;
2662 // for (unsigned int i = 0; i < Cb->size(0); ++i)
2663 // for (unsigned int j = i ; j < Cb->size(1); ++j)
2664 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*Cb)(i, j) - (*Ab)(i, j) / a1) < tol, true);
2665
2666 // Block = a * Simple
2667 *Cb = a * *A;
2668 for(unsigned int i = 0; i < Cb->size(0); ++i)
2669 for(unsigned int j = i ; j < Cb->size(1); ++j)
2670 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*Cb)(i, j) - a * (*A)(i, j)) < tol, true);
2671 *Cb = a1 * *A;
2672 for(unsigned int i = 0; i < Cb->size(0); ++i)
2673 for(unsigned int j = i ; j < Cb->size(1); ++j)
2674 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*Cb)(i, j) - a1 * (*A)(i, j)) < tol, true);
2675
2676 // *Cb = *A / a;
2677 // for (unsigned int i = 0; i < Cb->size(0); ++i)
2678 // for (unsigned int j = i ; j < Cb->size(1); ++j)
2679 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*Cb)(i, j) - (*A)(i, j) / a) < tol, true);
2680 // *Cb = *A / a1;
2681 // for (unsigned int i = 0; i < Cb->size(0); ++i)
2682 // for (unsigned int j = i ; j < Cb->size(1); ++j)
2683 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9: ", fabs((*Cb)(i, j) - (*A)(i, j) / a1) < tol, true);
2684 std::cout << "--> test operators9 ended with success." <<std::endl;
2685 }
2686
testOperators9Bis()2687 void SimpleMatrixTest::testOperators9Bis()
2688 {
2689 std::cout << "--> Test: operator9Bis." <<std::endl;
2690
2691 // C = a*A or A/a
2692
2693 double a = 2.2;
2694
2695 // Simple = a * Simple or Simple/a
2696 scal(a, *A, *C);
2697 for(unsigned int i = 0; i < C->size(0); ++i)
2698 for(unsigned int j = i ; j < C->size(1); ++j)
2699 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9Bis: ", fabs((*C)(i, j) - a * (*A)(i, j)) < tol, true);
2700
2701 scal(1.0 / a, *A, *C);
2702 for(unsigned int i = 0; i < C->size(0); ++i)
2703 for(unsigned int j = i ; j < C->size(1); ++j)
2704 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9Bis: ", fabs((*C)(i, j) - (*A)(i, j) / a) < tol, true);
2705 // Simple = a * Block
2706
2707 scal(a, *Ab, *C);
2708 for(unsigned int i = 0; i < C->size(0); ++i)
2709 for(unsigned int j = i ; j < C->size(1); ++j)
2710 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9Bis: ", fabs((*C)(i, j) - a * (*Ab)(i, j)) < tol, true);
2711
2712 scal(1.0 / a, *Ab, *C);
2713 for(unsigned int i = 0; i < C->size(0); ++i)
2714 for(unsigned int j = i ; j < C->size(1); ++j)
2715 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9Bis: ", fabs((*C)(i, j) - (*Ab)(i, j) / a) < tol, true);
2716
2717 // Block = a * Block
2718 scal(a, *Ab, *Cb);
2719 for(unsigned int i = 0; i < Cb->size(0); ++i)
2720 for(unsigned int j = i ; j < Cb->size(1); ++j)
2721 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9Bis: ", fabs((*Cb)(i, j) - a * (*Ab)(i, j)) < tol, true);
2722
2723 scal(1.0 / a, *Ab, *Cb);
2724 for(unsigned int i = 0; i < Cb->size(0); ++i)
2725 for(unsigned int j = i ; j < Cb->size(1); ++j)
2726 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9Bis: ", fabs((*Cb)(i, j) - (*Ab)(i, j) / a) < tol, true);
2727
2728 // Block = a * Simple
2729 scal(a, *A, *Cb);
2730 for(unsigned int i = 0; i < Cb->size(0); ++i)
2731 for(unsigned int j = i ; j < Cb->size(1); ++j)
2732 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9Bis: ", fabs((*Cb)(i, j) - a * (*A)(i, j)) < tol, true);
2733
2734 scal(1.0 / a, *A, *Cb);
2735 for(unsigned int i = 0; i < Cb->size(0); ++i)
2736 for(unsigned int j = i ; j < Cb->size(1); ++j)
2737 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9Bis: ", fabs((*Cb)(i, j) - (*A)(i, j) / a) < tol, true);
2738 std::cout << "--> test operators9Bis ended with success." <<std::endl;
2739 }
2740
testOperators9Ter()2741 void SimpleMatrixTest::testOperators9Ter()
2742 {
2743 std::cout << "--> Test: operator9Ter." <<std::endl;
2744
2745 // C += a*A or A/a
2746
2747 double a = 2.2;
2748 C->zero();
2749 // Simple = a * Simple or Simple/a
2750 scal(a, *A, *C, false);
2751 scal(a, *A, *C, false);
2752 for(unsigned int i = 0; i < C->size(0); ++i)
2753 for(unsigned int j = i ; j < C->size(1); ++j)
2754 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9Ter: ", fabs((*C)(i, j) - 2 * a * (*A)(i, j)) < tol, true);
2755
2756 // Simple = a * Block
2757 C->zero();
2758 scal(a, *Ab, *C, false);
2759 scal(a, *Ab, *C, false);
2760 for(unsigned int i = 0; i < C->size(0); ++i)
2761 for(unsigned int j = i ; j < C->size(1); ++j)
2762 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9Ter: ", fabs((*C)(i, j) - 2 * a * (*Ab)(i, j)) < tol, true);
2763
2764 // Block = a * Block
2765 Cb->zero();
2766 scal(a, *Ab, *Cb, false);
2767 scal(a, *Ab, *Cb, false);
2768 for(unsigned int i = 0; i < Cb->size(0); ++i)
2769 for(unsigned int j = i ; j < Cb->size(1); ++j)
2770 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9Ter: ", fabs((*Cb)(i, j) - 2 * a * (*Ab)(i, j)) < tol, true);
2771
2772 // Block = a * Simple
2773 Cb->zero();
2774 scal(a, *A, *Cb, false);
2775 scal(a, *A, *Cb, false);
2776 for(unsigned int i = 0; i < Cb->size(0); ++i)
2777 for(unsigned int j = i ; j < Cb->size(1); ++j)
2778 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators9Ter: ", fabs((*Cb)(i, j) - 2 * a * (*A)(i, j)) < tol, true);
2779
2780 std::cout << "--> test operators9Ter ended with success." <<std::endl;
2781 }
2782
testOperators10()2783 void SimpleMatrixTest::testOperators10()
2784 {
2785 std::cout << "--> Test: operator10." <<std::endl;
2786 double m = 2.2;
2787 int i = 3;
2788 SP::SiconosMatrix tmp1(new SimpleMatrix(*T));
2789 SP::SiconosMatrix res(new SimpleMatrix(3, 3, TRIANGULAR));
2790 *res = m ** tmp1;
2791 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getTriang() - tmp1->getTriang()*m) < tol, true);
2792 *res = i ** tmp1;
2793 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getTriang() - tmp1->getTriang()*i) < tol, true);
2794 // *res = *tmp1 / m;
2795 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getTriang() - tmp1->getTriang() / m) < tol, true);
2796 // *res = *tmp1 / i;
2797 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getTriang() - tmp1->getTriang() / i) < tol, true);
2798 std::cout << "--> test operators10 ended with success." <<std::endl;
2799 }
2800
testOperators11()2801 void SimpleMatrixTest::testOperators11()
2802 {
2803 std::cout << "--> Test: operator11." <<std::endl;
2804 double m = 2.2;
2805 int i = 3;
2806 SP::SiconosMatrix tmp1(new SimpleMatrix(*S));
2807 SP::SiconosMatrix res(new SimpleMatrix(3, 3, SYMMETRIC));
2808 *res = m ** tmp1;
2809 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getSym() - tmp1->getSym()*m) < tol, true);
2810 *res = i ** tmp1;
2811 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getSym() - tmp1->getSym()*i) < tol, true);
2812 // *res = *tmp1 / m;
2813 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getSym() - tmp1->getSym() / m) < tol, true);
2814 // *res = *tmp1 / i;
2815 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getSym() - tmp1->getSym() / i) < tol, true);
2816 std::cout << "--> test operator11 ended with success." <<std::endl;
2817 }
2818
testOperators12()2819 void SimpleMatrixTest::testOperators12()
2820 {
2821 std::cout << "--> Test: operator12." <<std::endl;
2822 double m = 2.2;
2823 int i = 3;
2824 SP::SiconosMatrix tmp1(new SimpleMatrix(*SP));
2825 SP::SiconosMatrix res(new SimpleMatrix(4, 4, Siconos::SPARSE));
2826 *res = m ** tmp1;
2827 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getSparse() - tmp1->getSparse()*m) < tol, true);
2828 *res = i ** tmp1;
2829 CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getSparse() - tmp1->getSparse()*i) < tol, true);
2830 // *res = *tmp1 / m;
2831 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getSparse() - tmp1->getSparse() / m) < tol, true);
2832 // *res = *tmp1 / i;
2833 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getSparse() - tmp1->getSparse() / i) < tol, true);
2834 std::cout << "--> test operators12 ended with success." <<std::endl;
2835 }
2836
testOperators13()2837 void SimpleMatrixTest::testOperators13()
2838 {
2839 std::cout << "--> Test: operator13." <<std::endl;
2840 // double m = 2.2;
2841 // int i = 3;
2842 // SP::SiconosMatrix tmp1(new SimpleMatrix(*Band);
2843 // SP::SiconosMatrix res(new SimpleMatrix(*Band);//4,4,BANDED,1,1);
2844 // *res = m * *tmp1;
2845 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getBanded()- tmp1->getBanded()*m)<tol, true);
2846 // *res = i ** tmp1;
2847 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getBanded()- tmp1->getBanded()*i)<tol, true);
2848 // *res = *tmp1 * m;
2849 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getBanded()- tmp1->getBanded()*m)<tol, true);
2850 // *res = *tmp1 * i;
2851 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getBanded()- tmp1->getBanded()*i)<tol, true);
2852 // *res = *tmp1 / m;
2853 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getBanded()- tmp1->getBanded()/m)<tol, true);
2854 // *res = *tmp1 / i;
2855 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testOperators: ", norm_inf(res->getBanded()- tmp1->getBanded()/i)<tol, true);
2856 std::cout << "--> test operators13 ended with success." <<std::endl;
2857 }
2858
testProd()2859 void SimpleMatrixTest::testProd() // y = A*x
2860 {
2861 std::cout << "--> Test: prod. mat-vect" <<std::endl;
2862
2863 SP::SiconosVector y(new SiconosVector(size));
2864 SP::SiconosVector x(new SiconosVector(size, 4.3));
2865 SP::SiconosVector x1(new SiconosVector(size - 2, 2.3));
2866 SP::SiconosVector x2(new SiconosVector(2, 3.1));
2867
2868 SP::BlockVector xB(new BlockVector(x1, x2));
2869 SP::BlockVector yB(new BlockVector(*xB));
2870 yB->zero();
2871
2872 // Matrix - vector product
2873
2874 // Simple = Simple * Simple
2875 *y = prod(*A, *x);
2876 double sum;
2877 for(unsigned int i = 0; i < size; ++i)
2878 {
2879 sum = 0;
2880 for(unsigned int j = 0; j < A->size(1); ++j)
2881 sum += (*A)(i, j) * (*x)(j);
2882 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd: ", fabs((*y)(i) - sum) < tol, true);
2883 }
2884 // Simple = Simple * Block
2885 // *y = prod(*A , *xB);
2886 // for (unsigned int i = 0; i< size; ++i)
2887 // {
2888 // sum = 0;
2889 // for (unsigned int j=0; j< A->size(1); ++j)
2890 // sum += (*A)(i,j)*(*xB)(j);
2891 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd: ", fabs((*y)(i) - sum)< tol, true);
2892 // }
2893
2894
2895 // Block = Simple * Simple
2896 *yB = prod(*A, *x);
2897 for(unsigned int i = 0; i < size; ++i)
2898 {
2899 sum = 0;
2900 for(unsigned int j = 0; j < A->size(1); ++j)
2901 sum += (*A)(i, j) * (*x)(j);
2902 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd: ", fabs((*yB)(i) - sum) < tol, true);
2903 }
2904
2905 // Block = Simple * Block
2906 // *yB = prod(*A ,*xB);
2907 // for (unsigned int i = 0; i< size; ++i)
2908 // {
2909 // sum = 0;
2910 // for (unsigned int j=0; j< A->size(1); ++j)
2911 // sum += (*A)(i,j)*(*xB)(j);
2912 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd: ", fabs((*yB)(i) - sum)< tol, true);
2913 // }
2914
2915 // Others or old stuff ...
2916
2917 SP::SiconosMatrix tmp2(new SimpleMatrix(*T));
2918 SP::SiconosMatrix tmp3(new SimpleMatrix(*S));
2919 SP::SiconosMatrix tmp4(new SimpleMatrix(*SP));
2920 SP::SiconosMatrix tmp5(new SimpleMatrix(*Band2));
2921 SP::SiconosVector v(new SiconosVector(3));
2922 (*v)(0) = 1;
2923 (*v)(1) = 2;
2924 (*v)(2) = 3;
2925 SP::SiconosVector vv(new SiconosVector(4));
2926 (*vv)(0) = 1;
2927 (*vv)(1) = 2;
2928 (*vv)(2) = 3;
2929 SparseVect * sv(new SparseVect(3));
2930 (*sv)(0) = 4;
2931 (*sv)(1) = 5;
2932 (*sv)(2) = 6;
2933 SparseVect * sv2(new SparseVect(4));
2934 (*sv2)(0) = 4;
2935 (*sv2)(1) = 5;
2936 (*sv2)(2) = 6;
2937 SP::SiconosVector w(new SiconosVector(*sv));
2938 SP::SiconosVector ww(new SiconosVector(*sv2));
2939 SP::SiconosVector res(new SiconosVector(4));
2940 SP::SiconosVector res2(new SiconosVector(3));
2941
2942 // Triang * ...
2943 *res2 = prod(*tmp2, *v);
2944 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd: ", norm_2(*res2->dense() - prod(tmp2->getTriang(), *v->dense())) < tol, true);
2945 *res2 = prod(*tmp2, *w);
2946 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd: ", norm_2(*res2->dense() - prod(tmp2->getTriang(), *w->sparse())) < tol, true);
2947 // Sym * ...
2948 *res2 = prod(*tmp3, *v);
2949 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd: ", norm_2(*res2->dense() - prod(tmp3->getSym(), *v->dense())) < tol, true);
2950 *res2 = prod(*tmp3, *w);
2951 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd: ", norm_2(*res2->dense() - prod(tmp3->getSym(), *w->sparse())) < tol, true);
2952 // Sparse * ...
2953 *res = prod(*tmp4, *vv);
2954 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd: ", norm_2(*res->dense() - prod(tmp4->getSparse(), *vv->dense())) < tol, true);
2955 *res = prod(*tmp4, *ww);
2956 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd: ", norm_2(*res->dense() - prod(tmp4->getSparse(), *ww->sparse())) < tol, true);
2957 // Triang * ...
2958 *res = prod(*tmp5, *v);
2959 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd: ", norm_2(*res->dense() - prod(tmp5->getBanded(), *v->dense())) < tol, true);
2960 *res = prod(*tmp5, *w);
2961 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd: ", norm_2(*res->dense() - prod(tmp5->getBanded(), *w->sparse())) < tol, true);
2962 std::cout << "--> test prod ended with success." <<std::endl;
2963
2964 delete sv2;
2965 delete sv;
2966 }
2967
testProdBis()2968 void SimpleMatrixTest::testProdBis()
2969 {
2970 std::cout << "--> Test: prod. mat-vect (bis)" <<std::endl;
2971
2972 SP::SiconosVector y(new SiconosVector(size));
2973 SP::SiconosVector x(new SiconosVector(size, 4.3));
2974 SP::SiconosVector x1(new SiconosVector(size - 2, 2.3));
2975 SP::SiconosVector x2(new SiconosVector(2, 3.1));
2976
2977 SP::BlockVector xB(new BlockVector(x1, x2));
2978 SP::BlockVector yB(new BlockVector(*xB));
2979 yB->zero();
2980
2981 // Matrix - vector product
2982
2983 // Simple = Simple * Simple
2984 prod(*A, *x, *y);
2985 double sum;
2986 for(unsigned int i = 0; i < size; ++i)
2987 {
2988 sum = 0;
2989 for(unsigned int j = 0; j < A->size(1); ++j)
2990 sum += (*A)(i, j) * (*x)(j);
2991 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdBis: ", fabs((*y)(i) - sum) < tol, true);
2992 }
2993 // Simple = Simple * Block
2994 prod(*A, *xB, *y);
2995 for(unsigned int i = 0; i < size; ++i)
2996 {
2997 sum = 0;
2998 for(unsigned int j = 0; j < A->size(1); ++j)
2999 sum += (*A)(i, j) * (*xB)(j);
3000 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdBis: ", fabs((*y)(i) - sum) < tol, true);
3001 }
3002
3003 // Block = Simple * Simple
3004 prod(*A, *x, *yB);
3005 for(unsigned int i = 0; i < size; ++i)
3006 {
3007 sum = 0;
3008 for(unsigned int j = 0; j < A->size(1); ++j)
3009 sum += (*A)(i, j) * (*x)(j);
3010 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdBis: ", fabs((*yB)(i) - sum) < tol, true);
3011 }
3012
3013 // Block = Simple * Block
3014 // prod(*A ,*xB,*yB);
3015 // for (unsigned int i = 0; i< size; ++i)
3016 // {
3017 // sum = 0;
3018 // for (unsigned int j=0; j< A->size(1); ++j)
3019 // sum += (*A)(i,j)*(*xB)(j);
3020 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdBis: ", fabs((*yB)(i) - sum)< tol, true);
3021 // }
3022 //
3023 // Others or old stuff ...
3024
3025 SP::SiconosMatrix tmp2(new SimpleMatrix(*T));
3026 SP::SiconosMatrix tmp3(new SimpleMatrix(*S));
3027 SP::SiconosMatrix tmp4(new SimpleMatrix(*SP));
3028 SP::SiconosMatrix tmp5(new SimpleMatrix(*Band2));
3029 SP::SiconosVector v(new SiconosVector(3));
3030 (*v)(0) = 1;
3031 (*v)(1) = 2;
3032 (*v)(2) = 3;
3033 SP::SiconosVector vv(new SiconosVector(4));
3034 (*vv)(0) = 1;
3035 (*vv)(1) = 2;
3036 (*vv)(2) = 3;
3037 SP::SparseVect sv(new SparseVect(3));
3038 (*sv)(0) = 4;
3039 (*sv)(1) = 5;
3040 (*sv)(2) = 6;
3041 SP::SparseVect sv2(new SparseVect(4));
3042 (*sv2)(0) = 4;
3043 (*sv2)(1) = 5;
3044 (*sv2)(2) = 6;
3045 SP::SiconosVector w(new SiconosVector(*sv));
3046 SP::SiconosVector ww(new SiconosVector(*sv2));
3047 SP::SiconosVector res(new SiconosVector(4));
3048 SP::SiconosVector res2(new SiconosVector(3));
3049
3050 // Triang * ...
3051 prod(*tmp2, *v, *res2);
3052 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdBis: ", norm_2(*res2->dense() - prod(tmp2->getTriang(), *v->dense())) < tol, true);
3053 prod(*tmp2, *w, *res2);
3054 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdBis: ", norm_2(*res2->dense() - prod(tmp2->getTriang(), *w->sparse())) < tol, true);
3055 // Sym * ...
3056 prod(*tmp3, *v, *res2);
3057 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdBis: ", norm_2(*res2->dense() - prod(tmp3->getSym(), *v->dense())) < tol, true);
3058 prod(*tmp3, *w, *res2);
3059 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdBis: ", norm_2(*res2->dense() - prod(tmp3->getSym(), *w->sparse())) < tol, true);
3060 // Sparse * ...
3061 prod(*tmp4, *vv, *res);
3062 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdBis: ", norm_2(*res->dense() - prod(tmp4->getSparse(), *vv->dense())) < tol, true);
3063 prod(*tmp4, *ww, *res);
3064 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdBis: ", norm_2(*res->dense() - prod(tmp4->getSparse(), *ww->sparse())) < tol, true);
3065 // Banded * ...
3066 prod(*tmp5, *v, *res);
3067 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdBis: ", norm_2(*res->dense() - prod(tmp5->getBanded(), *v->dense())) < tol, true);
3068 prod(*tmp5, *w, *res);
3069 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdBis: ", norm_2(*res->dense() - prod(tmp5->getBanded(), *w->sparse())) < tol, true);
3070 std::cout << "--> test prodBis ended with success." <<std::endl;
3071 }
3072
testProdTer()3073 void SimpleMatrixTest::testProdTer()
3074 {
3075 std::cout << "--> Test: prod. mat-vect (ter)" <<std::endl;
3076
3077 SP::SiconosVector y(new SiconosVector(size));
3078 SP::SiconosVector x(new SiconosVector(size, 4.3));
3079 SP::SiconosVector x1(new SiconosVector(size - 2, 2.3));
3080 SP::SiconosVector x2(new SiconosVector(2, 3.1));
3081
3082 SP::BlockVector xB(new BlockVector(x1, x2));
3083 SP::BlockVector yB(new BlockVector(*xB));
3084 yB->zero();
3085
3086 // Matrix - vector product
3087
3088 // Simple = Simple * Simple
3089 // axpy_prod(*A, *x, *y, true);
3090 // double sum;
3091 // for (unsigned int i = 0; i < size; ++i)
3092 // {
3093 // sum = 0;
3094 // for (unsigned int j = 0; j < A->size(1); ++j)
3095 // sum += (*A)(i, j) * (*x)(j);
3096 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", fabs((*y)(i) - sum) < tol, true);
3097 // }
3098
3099 //SP::SiconosVector backUp(new SiconosVector(*y));
3100 // // Simple += Simple * Simple
3101 // axpy_prod(*A, *x, *y, false);
3102 // for (unsigned int i = 0; i < size; ++i)
3103 // {
3104 // sum = 0;
3105 // for (unsigned int j = 0; j < A->size(1); ++j)
3106 // sum += (*A)(i, j) * (*x)(j);
3107 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", fabs((*y)(i) - sum - (*backUp)(i)) < tol, true);
3108 // }
3109
3110 // Simple = Simple * Block
3111 // axpy_prod(*A ,*xB,*y, true);
3112 // for (unsigned int i = 0; i< size; ++i)
3113 // {
3114 // sum = 0;
3115 // for (unsigned int j=0; j< A->size(1); ++j)
3116 // sum += (*A)(i,j)*(*xB)(j);
3117 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", fabs((*y)(i) - sum)< tol, true);
3118 // }
3119 //
3120 //*backUp = *y;
3121 // Simple += Simple * Block
3122 // axpy_prod(*A ,*xB,*y, false);
3123 // for (unsigned int i = 0; i< size; ++i)
3124 // {
3125 // sum = 0;
3126 // for (unsigned int j=0; j< A->size(1); ++j)
3127 // sum += (*A)(i,j)*(*xB)(j);
3128 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", fabs((*y)(i) - sum - (*backUp)(i))< tol, true);
3129 // }
3130 //
3131 // // Block = Simple * Simple
3132 // axpy_prod(*A ,*x,*yB, true);
3133 // for (unsigned int i = 0; i< size; ++i)
3134 // {
3135 // sum = 0;
3136 // for (unsigned int j=0; j< A->size(1); ++j)
3137 // sum += (*A)(i,j)*(*x)(j);
3138 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", fabs((*yB)(i) - sum)< tol, true);
3139 // }
3140 //
3141 // // Block += Simple * Simple
3142 // *backUp = *yB;
3143 // axpy_prod(*A ,*x,*yB, false);
3144 // for (unsigned int i = 0; i< size; ++i)
3145 // {
3146 // sum = 0;
3147 // for (unsigned int j=0; j< A->size(1); ++j)
3148 // sum += (*A)(i,j)*(*x)(j);
3149 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", fabs((*yB)(i) - sum - (*backUp)(i))< tol, true);
3150 // }
3151 //
3152 // // Block = Simple * Block
3153 // axpy_prod(*A ,*xB,*yB,true);
3154 // for (unsigned int i = 0; i< size; ++i)
3155 // {
3156 // sum = 0;
3157 // for (unsigned int j=0; j< A->size(1); ++j)
3158 // sum += (*A)(i,j)*(*xB)(j);
3159 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", fabs((*yB)(i) - sum)< tol, true);
3160 // }
3161 //
3162 // // Block += Simple * Block
3163 // *backUp = *yB;
3164 // axpy_prod(*A ,*xB,*yB,false);
3165 // for (unsigned int i = 0; i< size; ++i)
3166 // {
3167 // sum = 0;
3168 // for (unsigned int j=0; j< A->size(1); ++j)
3169 // sum += (*A)(i,j)*(*xB)(j);
3170 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", fabs((*yB)(i) - sum - (*backUp)(i))< tol, true);
3171 // }
3172 // // Others or old stuff ...
3173
3174 SP::SiconosMatrix tmp2(new SimpleMatrix(*T));
3175 SP::SiconosMatrix tmp3(new SimpleMatrix(*S));
3176 SP::SiconosMatrix tmp4(new SimpleMatrix(*SP));
3177 SP::SiconosMatrix tmp5(new SimpleMatrix(*Band2));
3178 SP::SiconosVector v(new SiconosVector(3));
3179 (*v)(0) = 1;
3180 (*v)(1) = 2;
3181 (*v)(2) = 3;
3182 SP::SiconosVector vv(new SiconosVector(4));
3183 (*vv)(0) = 1;
3184 (*vv)(1) = 2;
3185 (*vv)(2) = 3;
3186 SP::SparseVect sv(new SparseVect(3));
3187 (*sv)(0) = 4;
3188 (*sv)(1) = 5;
3189 (*sv)(2) = 6;
3190 SP::SparseVect sv2(new SparseVect(4));
3191 (*sv2)(0) = 4;
3192 (*sv2)(1) = 5;
3193 (*sv2)(2) = 6;
3194 SP::SiconosVector w(new SiconosVector(*sv));
3195 SP::SiconosVector ww(new SiconosVector(*sv2));
3196 SP::SiconosVector res(new SiconosVector(4));
3197 SP::SiconosVector res2(new SiconosVector(3));
3198
3199 // Triang * ...
3200 prod(*tmp2, *v, *res2);
3201 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", norm_2(*res2->dense() - prod(tmp2->getTriang(), *v->dense())) < tol, true);
3202 prod(*tmp2, *w, *res2);
3203 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", norm_2(*res2->dense() - prod(tmp2->getTriang(), *w->sparse())) < tol, true);
3204 // Sym * ...
3205 prod(*tmp3, *v, *res2);
3206 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", norm_2(*res2->dense() - prod(tmp3->getSym(), *v->dense())) < tol, true);
3207 prod(*tmp3, *w, *res2);
3208 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", norm_2(*res2->dense() - prod(tmp3->getSym(), *w->sparse())) < tol, true);
3209 // Sparse * ...
3210 prod(*tmp4, *vv, *res);
3211 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", norm_2(*res->dense() - prod(tmp4->getSparse(), *vv->dense())) < tol, true);
3212 prod(*tmp4, *ww, *res);
3213 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", norm_2(*res->dense() - prod(tmp4->getSparse(), *ww->sparse())) < tol, true);
3214 // Banded * ...
3215 prod(*tmp5, *v, *res);
3216 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", norm_2(*res->dense() - prod(tmp5->getBanded(), *v->dense())) < tol, true);
3217 prod(*tmp5, *w, *res);
3218 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProdTer: ", norm_2(*res->dense() - prod(tmp5->getBanded(), *w->sparse())) < tol, true);
3219 std::cout << "--> test prodTer ended with success." <<std::endl;
3220 }
3221
testProd4()3222 void SimpleMatrixTest::testProd4() // y += A*x
3223 {
3224 std::cout << "--> Test: prod. mat-vect (4)" <<std::endl;
3225
3226 SP::SiconosVector y(new SiconosVector(size));
3227 SP::SiconosVector x(new SiconosVector(size, 4.3));
3228 SP::SiconosVector x1(new SiconosVector(size - 2, 2.3));
3229 SP::SiconosVector x2(new SiconosVector(2, 3.1));
3230
3231 SP::BlockVector xB(new BlockVector(x1, x2));
3232 SP::BlockVector yB(new BlockVector(*xB));
3233 yB->zero();
3234
3235 // Matrix - vector product
3236
3237 // Simple = Simple * Simple
3238 y->zero();
3239 prod(*A, *x, *y, false);
3240 prod(*A, *x, *y, false);
3241 double sum;
3242 for(unsigned int i = 0; i < size; ++i)
3243 {
3244 sum = 0;
3245 for(unsigned int j = 0; j < A->size(1); ++j)
3246 sum += 2 * (*A)(i, j) * (*x)(j);
3247 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd4: ", fabs((*y)(i) - sum) < tol, true);
3248 }
3249 // Simple = Simple * Block
3250 y->zero();
3251 prod(*A, *xB, *y, false);
3252 prod(*A, *xB, *y, false);
3253 for(unsigned int i = 0; i < size; ++i)
3254 {
3255 sum = 0;
3256 for(unsigned int j = 0; j < A->size(1); ++j)
3257 sum += 2 * (*A)(i, j) * (*xB)(j);
3258 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd4: ", fabs((*y)(i) - sum) < tol, true);
3259 }
3260
3261 // Block = Simple * Simple
3262 yB->zero();
3263 prod(*A, *x, *yB, false);
3264 prod(*A, *x, *yB, false);
3265 for(unsigned int i = 0; i < size; ++i)
3266 {
3267 sum = 0;
3268 for(unsigned int j = 0; j < A->size(1); ++j)
3269 sum += 2 * (*A)(i, j) * (*x)(j);
3270 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd4: ", fabs((*yB)(i) - sum) < tol, true);
3271 }
3272
3273 // Block = Simple * Block
3274 yB->zero();
3275 // prod(*A ,*xB,*yB,false);
3276 // prod(*A ,*xB,*yB,false);
3277 // for (unsigned int i = 0; i< size; ++i)
3278 // {
3279 // sum = 0;
3280 // for (unsigned int j=0; j< A->size(1); ++j)
3281 // sum += 2*(*A)(i,j)*(*xB)(j);
3282 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd4: ", fabs((*yB)(i) - sum)< tol, true);
3283 // }
3284 std::cout << "--> test prod4 ended with success." <<std::endl;
3285 }
3286
testProd5()3287 void SimpleMatrixTest::testProd5() // y += a*A*x
3288 {
3289 std::cout << "--> Test: prod. mat-vect (5)" <<std::endl;
3290
3291 SP::SiconosVector y(new SiconosVector(size));
3292 SP::SiconosVector x(new SiconosVector(size, 4.3));
3293 SP::SiconosVector x1(new SiconosVector(size - 2, 2.3));
3294 SP::SiconosVector x2(new SiconosVector(2, 3.1));
3295
3296 SP::BlockVector xB(new BlockVector(x1, x2));
3297 SP::BlockVector yB(new BlockVector(*xB));
3298 yB->zero();
3299
3300 // Matrix - vector product
3301 double a = 3.0;
3302 // Simple = Simple * Simple
3303 y->zero();
3304 prod(a, *A, *x, *y, false);
3305 prod(a, *A, *x, *y, false);
3306 double sum;
3307 for(unsigned int i = 0; i < size; ++i)
3308 {
3309 sum = 0;
3310 for(unsigned int j = 0; j < A->size(1); ++j)
3311 sum += 2 * a * (*A)(i, j) * (*x)(j);
3312 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd5: ", fabs((*y)(i) - sum) < tol, true);
3313 }
3314 // Simple = Simple * Block
3315 y->zero();
3316 // prod(a,*A ,*xB,*y,false);
3317 // prod(a,*A ,*xB,*y,false);
3318 // for (unsigned int i = 0; i< size; ++i)
3319 // {
3320 // sum = 0;
3321 // for (unsigned int j=0; j< A->size(1); ++j)
3322 // sum += a*2*(*A)(i,j)*(*xB)(j);
3323 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd5: ", fabs((*y)(i) - sum)< tol, true);
3324 // }
3325
3326 // Block = Simple * Simple
3327 yB->zero();
3328 // prod(a,*A ,*x,*yB,false);
3329 // prod(a,*A ,*x,*yB,false);
3330 // for (unsigned int i = 0; i< size; ++i)
3331 // {
3332 // sum = 0;
3333 // for (unsigned int j=0; j< A->size(1); ++j)
3334 // sum += a*2*(*A)(i,j)*(*x)(j);
3335 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd5: ", fabs((*yB)(i) - sum)< tol, true);
3336 // }
3337 //
3338 // Block = Simple * Block
3339 yB->zero();
3340 // prod(a,*A ,*xB,*yB,false);
3341 // prod(a,*A ,*xB,*yB,false);
3342 // for (unsigned int i = 0; i< size; ++i)
3343 // {
3344 // sum = 0;
3345 // for (unsigned int j=0; j< A->size(1); ++j)
3346 // sum += a*2*(*A)(i,j)*(*xB)(j);
3347 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd5: ", fabs((*yB)(i) - sum)< tol, true);
3348 // }
3349 std::cout << "--> test prod5 ended with success." <<std::endl;
3350 }
3351
testProd6()3352 void SimpleMatrixTest::testProd6() // y += trans(A)*x
3353 {
3354 std::cout << "--> Test: prod. mat-vect (6)" <<std::endl;
3355
3356 SP::SiconosVector y(new SiconosVector(size));
3357 SP::SiconosVector x(new SiconosVector(size, 4.3));
3358 SP::SiconosVector x1(new SiconosVector(size - 2, 2.3));
3359 SP::SiconosVector x2(new SiconosVector(2, 3.1));
3360
3361 SP::BlockVector xB(new BlockVector(x1, x2));
3362 SP::BlockVector yB(new BlockVector(*xB));
3363 yB->zero();
3364
3365 SP::SiconosMatrix tmp(new SimpleMatrix(*A));
3366 tmp->trans();
3367 // Matrix - vector product
3368
3369 // Simple = Simple * Simple
3370 y->zero();
3371 prod(*x, *A, *y);
3372 prod(*x, *A, *y, false);
3373 double sum;
3374 for(unsigned int i = 0; i < size; ++i)
3375 {
3376 sum = 0;
3377 for(unsigned int j = 0; j < A->size(1); ++j)
3378 sum += 2 * (*tmp)(i, j) * (*x)(j);
3379 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd6: ", fabs((*y)(i) - sum) < tol, true);
3380 }
3381 // Simple = Simple * Block
3382 y->zero();
3383 // prod(*xB,*A,*y);
3384 // prod(*xB,*A,*y,false);
3385 //
3386 // for (unsigned int i = 0; i< size; ++i)
3387 // {
3388 // sum = 0;
3389 // for (unsigned int j=0; j< size; ++j)
3390 // sum += 2*(*tmp)(i,j)*(*xB)(j);
3391 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd6: ", fabs((*y)(i) - sum)< tol, true);
3392 // }
3393
3394 // Block = Simple * Simple
3395 yB->zero();
3396 prod(*x, *A, *yB);
3397 prod(*x, *A, *yB, false);
3398 for(unsigned int i = 0; i < size; ++i)
3399 {
3400 sum = 0;
3401 for(unsigned int j = 0; j < A->size(1); ++j)
3402 sum += 2 * (*tmp)(i, j) * (*x)(j);
3403 CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd6: ", fabs((*yB)(i) - sum) < tol, true);
3404 }
3405
3406 // Block = Simple * Block
3407 yB->zero();
3408 // prod(*xB,*A ,*yB);
3409 // prod(*xB,*A ,*yB,false);
3410 // for (unsigned int i = 0; i< size; ++i)
3411 // {
3412 // sum = 0;
3413 // for (unsigned int j=0; j< A->size(1); ++j)
3414 // sum += 2*(*tmp)(i,j)*(*xB)(j);
3415 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testProd6: ", fabs((*yB)(i) - sum)< tol, true);
3416 // }
3417 std::cout << "--> test prod6 ended with success." <<std::endl;
3418 }
3419
3420 // void SimpleMatrixTest::testGemv()
3421 // {
3422 // std::cout << "--> Test: gemv" <<std::endl;
3423
3424 // SP::SiconosVector y(new SiconosVector(size, 1.0));
3425 // SP::SiconosVector x(new SiconosVector(size, 4.3));
3426
3427 // SP::SiconosVector backUp(new SiconosVector(*y));
3428
3429 // double a = 2.3;
3430 // double b = 1.5;
3431 // double sum;
3432 // gemv(a, *A, *x, b, *y);
3433
3434 // for (unsigned int i = 0; i < size; ++i)
3435 // {
3436 // sum = b * (*backUp)(i);
3437 // for (unsigned int j = 0; j < A->size(1); ++j)
3438 // sum += a * (*A)(i, j) * (*x)(j) ;
3439 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testgemv: ", fabs((*y)(i) - sum) < tol, true);
3440 // }
3441
3442 // *y = *backUp;
3443 // gemvtranspose(a, *A, *x, b, *y);
3444 // for (unsigned int i = 0; i < size; ++i)
3445 // {
3446 // sum = b * (*backUp)(i);
3447 // for (unsigned int j = 0; j < A->size(0); ++j)
3448 // sum += a * (*A)(j, i) * (*x)(j);
3449 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testgemv (trans): ", fabs((*y)(i) - sum) < tol, true);
3450 // }
3451 // std::cout << "--> test gemv ended with success." <<std::endl;
3452 // }
3453
3454 // void SimpleMatrixTest::testGemm()
3455 // {
3456 // std::cout << "--> Test: gemm." <<std::endl;
3457
3458 // double a = 2.3;
3459 // double b = 1.5;
3460 // *C = *A;
3461 // SP::SiconosMatrix backUp(new SimpleMatrix(*C));
3462
3463 // gemm(a, *A, *B, b, *C);
3464 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testGemm: ", norm_inf(*C->dense() - a * prod(*A->dense(), *B->dense()) - b**backUp->dense()) < tol, true);
3465
3466 // *C = *backUp;
3467 // gemmtranspose(a, *A, *B, b, *C);
3468 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testGemm (trans): ", norm_inf(*C->dense() - a * prod(trans(*A->dense()), trans(*B->dense())) - b**backUp->dense()) < tol, true);
3469 // std::cout << "--> test gemm ended with success." <<std::endl;
3470 // }
3471
3472
testFromAndFillCSC()3473 void SimpleMatrixTest::testFromAndFillCSC()
3474 {
3475 std::cout << "Start SimpleMatrixTest::testFromAndFillCSC() "<< std::endl;
3476 SP::SiconosMatrix Sparse4(new SimpleMatrix(*SP4));
3477 Sparse4->updateNumericsMatrix();
3478 NumericsMatrix *NM = Sparse4->numericsMatrix();
3479 NM_display(NM);
3480 // NumericsMatrix *NM_1 = NM_create(4,4, NM_SPARSE);
3481
3482 SP::SiconosMatrix Sparse1(new SimpleMatrix(4,4,Siconos::SPARSE));
3483 Sparse1->fromCSC(NM_csc(NM));
3484 Sparse1->displayExpert();
3485
3486 NumericsMatrix *NM_1 = NM_create(NM_SPARSE, 4,4);
3487 NM_1->matrix2->origin = NSM_CSC;
3488 NM_csc_alloc(NM_1, Sparse4->nnz());
3489 Sparse4->fillCSC(NM_csc(NM_1));
3490 //NM_display(NM_1); --> Note FP : fails when exiting the function ... To be investigating ...
3491 NM_1 = NM_free(NM_1);
3492 std::cout << "End SimpleMatrixTest::testFromAndFillCSC() "<< std::endl;
3493
3494 }
testPLUFactorizationInPlace()3495 void SimpleMatrixTest::testPLUFactorizationInPlace()
3496 {
3497 std::cout << "--> Test: PLUFactorizationInPlace." <<std::endl;
3498
3499 SP::SiconosMatrix Dense(new SimpleMatrix(*D));
3500 Dense->display();
3501 Dense->PLUFactorizationInPlace();
3502 Dense->display();
3503 //CPPUNIT_ASSERT_EQUAL_MESSAGE("testPLUFactorizationInPlace: ", < tol, true);
3504
3505 SP::SiconosMatrix Sparse(new SimpleMatrix(4,4,SPARSE));
3506 Sparse->eye();
3507 Sparse->display();
3508 Sparse->PLUFactorizationInPlace();
3509 Sparse->display();
3510
3511 SP::SimpleMatrix Sparse2(new SimpleMatrix(*SP4));
3512 DEBUG_EXPR(Sparse2->display(););
3513 Sparse2->PLUFactorizationInPlace();
3514 DEBUG_EXPR(Sparse2->display(););
3515
3516 std::cout << "--> test PLUFactorizationInPlace ended with success." <<std::endl;
3517 }
3518
testFactorize()3519 void SimpleMatrixTest::testFactorize()
3520 {
3521 std::cout << "--> Test: Factorize (LU)." <<std::endl;
3522
3523 SP::SiconosMatrix Dense(new SimpleMatrix(*D));
3524 Dense->display();
3525 Dense->Factorize();
3526 Dense->display();
3527 //CPPUNIT_ASSERT_EQUAL_MESSAGE("testFactorize: ", < tol, true);
3528
3529
3530 SP::SimpleMatrix Sparse(new SimpleMatrix(4,4,SPARSE));
3531 Sparse->eye();
3532 Sparse->display();
3533 Sparse->Factorize();
3534
3535
3536 Sparse.reset(new SimpleMatrix(*SP3));
3537 Sparse->display();
3538 Sparse->displayExpert();
3539 Sparse->Factorize();
3540 Sparse->display();
3541
3542
3543
3544 // Other types than SPARSE are not working since fillCSC is not implemented.
3545 // std::cout << "--> Test: Factorize -- Triangle" <<std::endl;
3546 // SP::SimpleMatrix Triangle(new SimpleMatrix(*T2));
3547 // Triangle->display();
3548 // Triangle->displayExpert();
3549 // Triangle->Factorize();
3550 // Triangle->display();
3551
3552
3553
3554 /* A last column full of zero caused memory corruption in cs_lusol
3555 since ublas does not fill the last entry correctly*/
3556 // Sparse.reset(new SimpleMatrix(*SP2));
3557 // Sparse->display();
3558 // Sparse->displayExpert();
3559 // Sparse->PLUFactorizationInPlace();
3560 // Sparse->display();
3561
3562 std::cout << "--> Test: Factorize (Cholesky)." <<std::endl;
3563
3564 Dense.reset(new SimpleMatrix(*D));
3565 /* conpute DD^T */
3566 SP::SiconosMatrix DenseT(new SimpleMatrix(*Dense));
3567 DenseT->trans();
3568 SP::SiconosMatrix DDT(new SimpleMatrix(Dense->size(0), Dense->size(1)));
3569 prod(*Dense, *DenseT, *DDT);
3570 DDT->display();
3571 DDT->setIsSymmetric(true);
3572 DDT->setIsPositiveDefinite(true);
3573 DDT->Factorize();
3574
3575 DDT->display();
3576 //CPPUNIT_ASSERT_EQUAL_MESSAGE("testFactorize: ", < tol, true);
3577
3578 std::cout << "--> test Factorize ended with success." <<std::endl;
3579 }
testSolve()3580 void SimpleMatrixTest::testSolve()
3581 {
3582 std::cout << "\n--> Test: Solve. Dense. LU." <<std::endl;
3583
3584 // Test dense matrix
3585 SP::SiconosMatrix Dense(new SimpleMatrix(*D));
3586 SP::SiconosVector b (new SiconosVector(Dense->size(0)));
3587 for( int i =0; i <Dense->size(0); i++)
3588 {
3589 (*b)(i)=1.0;
3590 }
3591 SP::SiconosVector backup (new SiconosVector(*b));
3592 SP::SiconosMatrix D_backup (new SimpleMatrix(*Dense));
3593 Dense->display();
3594 Dense->Solve(*b);
3595 Dense->display();
3596 b->display();
3597 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*D_backup,*b) - *backup).norm2() < tol, true);
3598
3599 // // Test dense matrix and sparse rhs
3600 // Dense.reset(new SimpleMatrix(*D));
3601 // SP::SiconosVector b_sparse (new SiconosVector(Dense->size(0), SPARSE));
3602 // for( int i =0; i <Dense->size(0); i++)
3603 // {
3604 // (*b_sparse)(i)=1.0;
3605 // }
3606 // backup.reset(new SiconosVector(*b_sparse));
3607 // Dense->Solve(*b_sparse);
3608 // b_sparse->display();
3609
3610 // CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*D_backup,*b_sparse) - *backup).norm2() < tol, true);
3611
3612
3613 std::cout << "\n\n--> Test: Solve. Dense. Cholesky." <<std::endl;
3614 Dense.reset(new SimpleMatrix(*D));
3615 b.reset(new SiconosVector(Dense->size(0)));
3616 for( int i =0; i <Dense->size(0); i++)
3617 {
3618 (*b)(i)=1.0;
3619 }
3620 SP::SiconosMatrix DenseT(new SimpleMatrix(*Dense));
3621 DenseT->trans();
3622 SP::SiconosMatrix DDT(new SimpleMatrix(Dense->size(0), Dense->size(1)));
3623 prod(*Dense, *DenseT, *DDT);
3624 SP::SiconosMatrix DDT_backup (new SimpleMatrix(*DDT));
3625 DDT->setIsSymmetric(true);
3626 DDT->setIsPositiveDefinite(true);
3627 DDT->Solve(*b);
3628 b->display();
3629 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*DDT_backup,*b) - *backup).norm2() < tol, true);
3630
3631
3632 std::cout << "\n\n--> Test: Solve. Sparse. LU." <<std::endl;
3633
3634 // Test sparse matrix identity
3635 SP::SimpleMatrix Sparse(new SimpleMatrix(4,4,SPARSE));
3636 SP::SimpleMatrix Sparse_backup(new SimpleMatrix(4,4,SPARSE));
3637 b.reset(new SiconosVector(Sparse->size(0)));
3638 for( int i =0; i <Sparse->size(0); i++)
3639 {
3640 (*b)(i)=1.0;
3641 }
3642 backup.reset (new SiconosVector(*b));
3643 Sparse_backup->eye();
3644 Sparse->eye();
3645 Sparse->Solve(*b);
3646 b->display();
3647 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*Sparse_backup,*b) - *backup).norm2() < tol, true);
3648
3649 // test sparse matrix 3x3
3650 Sparse.reset(new SimpleMatrix(*SP3));
3651 b.reset(new SiconosVector(Sparse->size(0)));
3652 for( int i =0; i <Sparse->size(0); i++)
3653 {
3654 (*b)(i)=1.0;
3655 }
3656 backup.reset (new SiconosVector(*b));
3657 Sparse->Solve(*b);
3658 b->display();
3659 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*Sparse,*b) - *backup).norm2() < tol, true);
3660
3661 // Solve again with another r.h.s.
3662 b.reset(new SiconosVector(Sparse->size(0)));
3663 for( int i =0; i <Sparse->size(0); i++)
3664 {
3665 (*b)(i)=2.0;
3666 }
3667 backup.reset (new SiconosVector(*b));
3668 Sparse->Solve(*b);
3669 b->display();
3670 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*Sparse,*b) - *backup).norm2() < tol, true);
3671
3672
3673 // test sparse matrix 4x4 SP4
3674 Sparse.reset(new SimpleMatrix(*SP4));
3675 b.reset(new SiconosVector(Sparse->size(0)));
3676 for( int i =0; i <Sparse->size(0); i++)
3677 {
3678 (*b)(i)=1.0;
3679 }
3680 backup.reset (new SiconosVector(*b));
3681 Sparse->Solve(*b);
3682 b->display();
3683 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*Sparse,*b) - *backup).norm2() < tol, true);
3684
3685
3686 std::cout << "\n\n--> Test: Solve. Sparse. LU. Sparse rhs" <<std::endl;
3687
3688
3689 // test sparse matrix 3x3 sparse RhS. trivial solution (Id)
3690 Sparse.reset(new SimpleMatrix(*SP3));
3691 SP::SimpleMatrix Sparse_rhs (new SimpleMatrix(*SP3));
3692 Sparse->Solve(*Sparse_rhs);
3693 // std::cout << "Sparse_rhs :" << std::endl;
3694 // Sparse_rhs->display();
3695 // std::cout << "Sparse :" << std::endl;
3696 // Sparse->display();
3697 // std::cout << "A A^{-1}" << std::endl;
3698 // (prod(*Sparse,*Sparse_rhs)).display();
3699 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*Sparse,*Sparse_rhs) - *Sparse ).normInf() < tol, true);
3700
3701 // test sparse matrix 3x3 sparse RhS. inverse
3702 Sparse.reset(new SimpleMatrix(*SP3));
3703 Sparse_rhs.reset(new SimpleMatrix(3,3));
3704 Sparse_rhs->eye();
3705 Sparse->Solve(*Sparse_rhs);
3706 SP::SiconosMatrix Id (new SimpleMatrix(3,3));
3707 Id->eye();
3708
3709 // Sparse_rhs->display();
3710 // std::cout << "A A^{-1}" << std::endl;
3711
3712 // (prod(*Sparse,*Sparse_rhs)).display();
3713 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*Sparse,*Sparse_rhs) - *Id ).normInf() < tol, true);
3714
3715
3716
3717 std::cout << "\n\n--> Test: Solve. Sparse. Cholesky." <<std::endl;
3718
3719 // Test sparse matrix identity
3720 Sparse.reset(new SimpleMatrix(4,4,SPARSE));
3721 Sparse_backup.reset(new SimpleMatrix(4,4,SPARSE));
3722 b.reset(new SiconosVector(Sparse->size(0)));
3723 for( int i =0; i <Sparse->size(0); i++)
3724 {
3725 (*b)(i)=1.0;
3726 }
3727 backup.reset (new SiconosVector(*b));
3728 Sparse_backup->eye();
3729 Sparse->eye();
3730 Sparse->setIsSymmetric(true);
3731 Sparse->setIsPositiveDefinite(true);
3732 Sparse->Solve(*b);
3733 b->display();
3734 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*Sparse_backup,*b) - *backup).norm2() < tol, true);
3735
3736 // test sparse matrix 3x3
3737 Sparse.reset(new SimpleMatrix(*SP3));
3738 b.reset(new SiconosVector(Sparse->size(0)));
3739 for( int i =0; i <Sparse->size(0); i++)
3740 {
3741 (*b)(i)=1.0;
3742 }
3743 backup.reset (new SiconosVector(*b));
3744 SP::SimpleMatrix SparseT (new SimpleMatrix(*SP3));
3745 SparseT->trans();
3746 SP::SimpleMatrix SST (new SimpleMatrix(Sparse->size(0), Sparse->size(1), SPARSE));
3747 prod(*Sparse, *SparseT, *SST);
3748 SST->setIsSymmetric(true);
3749 SST->setIsPositiveDefinite(true);
3750 SST->Solve(*b);
3751 b->display();
3752 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*SST,*b) - *backup).norm2() < tol, true);
3753
3754 // Solve again with another r.h.s.
3755 b.reset(new SiconosVector(Sparse->size(0)));
3756 for( int i =0; i <Sparse->size(0); i++)
3757 {
3758 (*b)(i)=2.0;
3759 }
3760 backup.reset (new SiconosVector(*b));
3761 SST->Solve(*b);
3762 b->display();
3763 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*SST,*b) - *backup).norm2() < tol, true);
3764
3765 // test sparse matrix 4x4 SP4
3766 Sparse.reset(new SimpleMatrix(*SP4));
3767 b.reset(new SiconosVector(Sparse->size(0)));
3768 for( int i =0; i <Sparse->size(0); i++)
3769 {
3770 (*b)(i)=1.0;
3771 }
3772 backup.reset (new SiconosVector(*b));
3773 SparseT.reset(new SimpleMatrix(*SP4));
3774 SparseT->trans();
3775 SST.reset (new SimpleMatrix(Sparse->size(0), Sparse->size(1), SPARSE));
3776 prod(*Sparse, *SparseT, *SST);
3777 SST->setIsSymmetric(true);
3778 SST->setIsPositiveDefinite(true);
3779 SST->Solve(*b);
3780 b->display();
3781 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*SST,*b) - *backup).norm2() < tol, true);
3782
3783
3784 std::cout << "\n\n--> Test: Solve. Sparse. Cholesky. Sparse rhs" <<std::endl;
3785
3786 // test sparse matrix 3x3 sparse RhS. trivial solution (Id)
3787 Sparse.reset(new SimpleMatrix(*SP3));
3788 SparseT.reset(new SimpleMatrix(*SP3));
3789 SparseT->trans();
3790 SST.reset (new SimpleMatrix(Sparse->size(0), Sparse->size(1), SPARSE));
3791 prod(*Sparse, *SparseT, *SST);
3792 Sparse_rhs.reset (new SimpleMatrix(*SST));
3793 // std::cout << "SST" << std::endl;
3794 // SST->display();
3795 SST->setIsSymmetric(true);
3796 SST->setIsPositiveDefinite(true);
3797 SST->Solve(*Sparse_rhs);
3798 // Sparse_rhs->display();
3799 // std::cout << "A A^{-1}" << std::endl;
3800 // (prod(*SST,*Sparse_rhs)).display();
3801 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*SST,*Sparse_rhs) - *SST ).normInf() < tol, true);
3802
3803 // test sparse matrix 3x3 sparse RhS. inverse
3804 SST.reset (new SimpleMatrix(Sparse->size(0), Sparse->size(1), SPARSE));
3805 prod(*Sparse, *SparseT, *SST);
3806 Sparse_rhs.reset(new SimpleMatrix(3,3, SPARSE));
3807 Sparse_rhs->eye();
3808 // std::cout << "SST" << std::endl;
3809 // SST->display();
3810 SST->setIsSymmetric(true);
3811 SST->setIsPositiveDefinite(true);
3812 SST->Solve(*Sparse_rhs);
3813 // Sparse_rhs->display();
3814
3815 // Sparse_rhs->display();
3816 // std::cout << "A A^{-1}" << std::endl;
3817
3818 // (prod(*SST,*Sparse_rhs)).display();
3819 CPPUNIT_ASSERT_EQUAL_MESSAGE("testSolve: ", (prod(*SST,*Sparse_rhs) - *Id ).normInf() < tol, true);
3820
3821
3822
3823 std::cout << "--> test Solve ended with success." <<std::endl;
3824 }
3825
3826
3827
End()3828 void SimpleMatrixTest::End()
3829 {
3830 std::cout << "======================================" <<std::endl;
3831 std::cout << " ===== End of SimpleMatrix Tests ===== " <<std::endl;
3832 std::cout << "======================================" <<std::endl;
3833 }
3834