1 /*
2 * Tiny Vector Matrix Library
3 * Dense Vector Matrix Libary of Tiny size using Expression Templates
4 *
5 * Copyright (C) 2001 - 2006 Olaf Petzold <opetzold@users.sourceforge.net>
6 *
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 *
21 * $Id: TestMathOps.h,v 1.2 2006-11-20 20:22:18 opetzold Exp $
22 */
23
24 #ifndef TVMET_TEST_MATHOPS_H
25 #define TVMET_TEST_MATHOPS_H
26
27 #include <cppunit/extensions/HelperMacros.h>
28
29 #include <tvmet/Vector.h>
30 #include <tvmet/Matrix.h>
31 #include <tvmet/util/General.h>
32
33 template <class T>
34 class TestMathOps : public CppUnit::TestFixture
35 {
36 CPPUNIT_TEST_SUITE( TestMathOps );
37 CPPUNIT_TEST( ScalarAssign );
38 CPPUNIT_TEST( Assign );
39 CPPUNIT_TEST( ScalarOps );
40 CPPUNIT_TEST( Ops1 );
41 CPPUNIT_TEST( Ops2 );
42 CPPUNIT_TEST( VectorOps );
43 CPPUNIT_TEST( VectorOps2 );
44 CPPUNIT_TEST( VectorNorm2 );
45 CPPUNIT_TEST( MatrixOps );
46 CPPUNIT_TEST( MatrixVector1 );
47 CPPUNIT_TEST( MatrixVector2 );
48 CPPUNIT_TEST( MatrixTransMatrix );
49 CPPUNIT_TEST( MatrixTransVector );
50 CPPUNIT_TEST( MatrixRowVector );
51 CPPUNIT_TEST( MatrixColVector );
52 CPPUNIT_TEST( MatrixDiagVector );
53 CPPUNIT_TEST( MatrixMatrixVector );
54 CPPUNIT_TEST_SUITE_END();
55
56 private:
57 typedef tvmet::Vector<T, 3> vector_type;
58 typedef tvmet::Matrix<T, 3, 3> matrix_type;
59
60 public:
TestMathOps()61 TestMathOps()
62 : vZero(0), vOne(1), mZero(0), mOne(1), scalar(10) { }
63
64 public: // cppunit interface
65 /** cppunit hook for fixture set up. */
66 void setUp();
67
68 /** cppunit hook for fixture tear down. */
69 void tearDown();
70
71 protected:
72 void ScalarAssign();
73 void Assign();
74 void ScalarOps();
75 void Ops1();
76 void Ops2();
77 void VectorOps();
78 void VectorOps2();
79 void VectorNorm2();
80 void MatrixOps();
81 void MatrixVector1();
82 void MatrixVector2();
83 void MatrixTransMatrix();
84 void MatrixTransVector();
85 void MatrixRowVector();
86 void MatrixColVector();
87 void MatrixDiagVector();
88 void MatrixMatrixVector();
89
90 private:
91 const vector_type vZero;
92 const vector_type vOne;
93 vector_type v1, v1b;
94 vector_type vBig; /**< vector 10x bigger than v1 */
95
96 private:
97 const matrix_type mZero;
98 const matrix_type mOne;
99 matrix_type m1, m1b;
100 matrix_type mBig; /**< matrix 10x bigger than m1 */
101
102 private:
103 vector_type m1_r0, m1_r1, m1_r2; // row vectors
104 vector_type m1_c0, m1_c1, m1_c2; // col vectors
105
106 private:
107 const T scalar;
108 };
109
110 /*****************************************************************************
111 * Implementation
112 ****************************************************************************/
113
114 /*
115 * cppunit part
116 */
117 template <class T>
setUp()118 void TestMathOps<T>::setUp() {
119 v1 = 1,2,3;
120 v1b = v1; // same as v1, cctor test done in checkInternal
121 vBig = 10,20,30;
122
123 m1 = 1,4,7,
124 2,5,8,
125 3,6,9;
126
127 m1_r0 = 1,4,7;
128 m1_r1 = 2,5,8;
129 m1_r2 = 3,6,9;
130
131 m1_c0 = 1,2,3;
132 m1_c1 = 4,5,6;
133 m1_c2 = 7,8,9;
134
135 m1b = m1; // same as m1, cctor test done in checkInternal
136
137 mBig = 10,40,70,
138 20,50,80,
139 30,60,90;
140
141 }
142
143 template <class T>
tearDown()144 void TestMathOps<T>::tearDown() {
145
146 }
147
148 /*
149 * regressions
150 */
151 template <class T>
152 void
ScalarAssign()153 TestMathOps<T>::ScalarAssign() {
154 {
155 vector_type t1(v1), t2(v1), t3(v1), t4(vBig);
156
157 t1 += scalar;
158 t2 -= scalar;
159 t3 *= scalar;
160 t4 /= scalar;
161
162 CPPUNIT_ASSERT(t1(0) == (v1(0)+scalar) && t1(1) == (v1(1)+scalar) && t1(2) == (v1(2)+scalar));
163 CPPUNIT_ASSERT(t2(0) == (v1(0)-scalar) && t2(1) == (v1(1)-scalar) && t2(2) == (v1(2)-scalar));
164 CPPUNIT_ASSERT( all_elements(t3 == vBig) );
165 CPPUNIT_ASSERT( all_elements(t4 == v1) );
166 }
167 {
168 matrix_type t1(m1), t2(m1), t3(m1), t4(mBig);
169
170 t1 += scalar;
171 t2 -= scalar;
172 t3 *= scalar;
173 t4 /= scalar;
174
175 CPPUNIT_ASSERT(t1(0,0) == (m1(0,0)+scalar) && t1(0,1) == (m1(0,1)+scalar) && t1(0,2) == (m1(0,2)+scalar) &&
176 t1(1,0) == (m1(1,0)+scalar) && t1(1,1) == (m1(1,1)+scalar) && t1(1,2) == (m1(1,2)+scalar) &&
177 t1(2,0) == (m1(2,0)+scalar) && t1(2,1) == (m1(2,1)+scalar) && t1(2,2) == (m1(2,2)+scalar));
178 CPPUNIT_ASSERT(t2(0,0) == (m1(0,0)-scalar) && t2(0,1) == (m1(0,1)-scalar) && t2(0,2) == (m1(0,2)-scalar) &&
179 t2(1,0) == (m1(1,0)-scalar) && t2(1,1) == (m1(1,1)-scalar) && t2(1,2) == (m1(1,2)-scalar) &&
180 t2(2,0) == (m1(2,0)-scalar) && t2(2,1) == (m1(2,1)-scalar) && t2(2,2) == (m1(2,2)-scalar));
181 CPPUNIT_ASSERT( all_elements(t3 == mBig) );
182 CPPUNIT_ASSERT( all_elements(t4 == m1) );
183 }
184 }
185
186 template <class T>
187 void
Assign()188 TestMathOps<T>::Assign() {
189 {
190 vector_type t1(vZero), t2(v1), t3(v1);
191
192 t1 += v1;
193 t2 -= v1;
194 t3 *= v1;
195
196 CPPUNIT_ASSERT( all_elements(t1 == v1) );
197 CPPUNIT_ASSERT( all_elements(t2 == vZero) );
198 CPPUNIT_ASSERT(t3(0) == (v1(0)*v1(0)) && t3(1) == (v1(1)*v1(1)) && t3(2) == (v1(2)*v1(2)));
199 }
200 {
201 matrix_type t1(mZero), t2(m1), t3(m1);
202
203 t1 += m1;
204 t2 -= m1;
205
206 CPPUNIT_ASSERT( all_elements(t1 == m1) );
207 CPPUNIT_ASSERT( all_elements(t2 == mZero) );
208 }
209 }
210
211 template <class T>
212 void
ScalarOps()213 TestMathOps<T>::ScalarOps() {
214 {
215 vector_type t1(v1), t2(v1), t3(v1), t4(vBig);
216 vector_type r1(v1), r2(v1);
217 r1 += scalar;
218 r2 -= scalar;
219
220 t1 = t1 + scalar;
221 t2 = t2 - scalar;
222 t3 = t3 * scalar;
223 t4 = t4 / scalar;
224
225 CPPUNIT_ASSERT( all_elements(t1 == r1) );
226 CPPUNIT_ASSERT( all_elements(t2 == r2) );
227 CPPUNIT_ASSERT( all_elements(t3 == vBig) );
228 CPPUNIT_ASSERT( all_elements(t4 == v1) );
229 }
230 {
231 matrix_type t1(m1), t2(m1), t3(m1), t4(mBig);
232 matrix_type r1(m1), r2(m1);
233 r1 += scalar;
234 r2 -= scalar;
235
236 t1 = t1 + scalar;
237 t2 = t2 - scalar;
238 t3 = t3 * scalar;
239 t4 = t4 / scalar;
240
241 CPPUNIT_ASSERT( all_elements(t1 == r1) );
242 CPPUNIT_ASSERT( all_elements(t2 == r2) );
243 CPPUNIT_ASSERT( all_elements(t3 == mBig) );
244 CPPUNIT_ASSERT( all_elements(t4 == m1) );
245 }
246 }
247
248 template <class T>
249 void
Ops1()250 TestMathOps<T>::Ops1() {
251 {
252 vector_type t1(0), t2(0), t3(0);
253 vector_type r(v1);
254 r *= v1;
255
256 t1 = v1 + v1;
257 t2 = v1 - v1;
258 t3 = v1 * v1;
259
260 CPPUNIT_ASSERT( all_elements(t1 == T(2)*v1) );
261 CPPUNIT_ASSERT( all_elements(t2 == vZero) );
262 CPPUNIT_ASSERT( all_elements(t3 == r) );
263 }
264 {
265 matrix_type t1(0), t2(0);
266 t1 = m1 + m1;
267 t2 = m1 - m1;
268
269 CPPUNIT_ASSERT( all_elements(t1 == T(2)*m1) );
270 CPPUNIT_ASSERT( all_elements(t2 == mZero) );
271 }
272 }
273
274 template <class T>
275 void
Ops2()276 TestMathOps<T>::Ops2() {
277 const vector_type vMinusOne(-1);
278 const matrix_type mMinusOne(-1);
279
280 // negate operator
281 {
282 vector_type t1, t2;
283
284 t1 = abs(v1);
285 CPPUNIT_ASSERT( all_elements(t1 == v1) );
286
287 t1 = -vOne;
288 CPPUNIT_ASSERT( all_elements(t1 == vMinusOne) );
289 }
290 {
291 matrix_type t1, t2;
292
293 t1 = abs(m1);
294 CPPUNIT_ASSERT( all_elements(t1 == m1) );
295
296 t1 = -mOne;
297 CPPUNIT_ASSERT( all_elements(t1 == mMinusOne) );
298
299 }
300 }
301
302 template <class T>
303 void
VectorOps()304 TestMathOps<T>::VectorOps() {
305
306 }
307
308 template <class T>
309 void
VectorOps2()310 TestMathOps<T>::VectorOps2() {
311 }
312
313 template <class T>
314 void
VectorNorm2()315 TestMathOps<T>::VectorNorm2() {
316 // casts for int vectors, as well as for complex<> since
317 // norm2 returns sum_type
318 CPPUNIT_ASSERT( norm2(v1) == static_cast<T>(std::sqrt(14.0)));
319 }
320
321 template <class T>
322 void
MatrixOps()323 TestMathOps<T>::MatrixOps() {
324 matrix_type t1, t2, t3;
325 matrix_type r1, r2, r3;
326
327 tvmet::util::Gemm(m1, m1, r1);
328 tvmet::util::Gemm(m1, mBig, r2);
329 tvmet::util::Gemm(mBig, m1, r3);
330 CPPUNIT_ASSERT( all_elements(r2 == r3) );
331
332 t1 = m1 * m1;
333 CPPUNIT_ASSERT( all_elements(t1 == r1) );
334
335 t2 = m1 * mBig;
336 CPPUNIT_ASSERT( all_elements(t2 == r2) );
337
338 t3 = mBig * m1;
339 CPPUNIT_ASSERT( all_elements(t3 == r3) );
340
341 t3 = trans(t1);
342 CPPUNIT_ASSERT( any_elements(t3 != t1) ); // XXX very simple test
343 t2 = trans(t3);
344 CPPUNIT_ASSERT( all_elements(t1 == t2) );
345
346 // trace return sum_type, therefore the cast for complex<>
347 CPPUNIT_ASSERT( static_cast<T>(trace(m1)) == static_cast<T>(15) );
348 }
349
350 template <class T>
351 void
MatrixVector1()352 TestMathOps<T>::MatrixVector1() {
353
354 vector_type t1, t2;
355 vector_type vr1(0), vr2(0); // clear it before use due to util::Gemv algo
356
357 // Matrix-Vector
358 tvmet::util::Gemv(m1, v1, vr1);
359 tvmet::util::Gemv(mBig, vBig, vr2);
360
361 t1 = m1 * v1;
362 t2 = mBig * vBig;
363
364 CPPUNIT_ASSERT( all_elements(t1 == vr1) );
365 CPPUNIT_ASSERT( all_elements(t2 == vr2) );
366 }
367
368 template <class T>
369 void
MatrixVector2()370 TestMathOps<T>::MatrixVector2() {
371
372 vector_type t1, t2;
373 vector_type vr(0), v2(0); // clear it before use due to util::Gemv algo
374
375 // Matrix-XprVector
376 v2 = v1 * vBig;
377 tvmet::util::Gemv(m1, v2, vr);
378
379 t1 = m1 * (v1*vBig);
380
381 CPPUNIT_ASSERT( all_elements(t1 == vr) );
382 }
383
384 template <class T>
385 void
MatrixTransMatrix()386 TestMathOps<T>::MatrixTransMatrix() {
387 // greatings to
388 {
389 matrix_type m1t, Mr, M2;
390
391 // trans() and prod() is checked before!
392 m1t = trans(m1);
393 Mr = prod(m1t, mBig);
394
395 M2 = MtM_prod(m1, mBig);
396
397 CPPUNIT_ASSERT( all_elements(Mr == M2) );
398 }
399 }
400
401 template <class T>
402 void
MatrixTransVector()403 TestMathOps<T>::MatrixTransVector() {
404 // greatings to
405 {
406 matrix_type Mt;
407 vector_type vr, y;
408
409 // trans() and prod() is checked before!
410 Mt = trans(m1);
411 vr = Mt*v1;
412 y = Mtx_prod(m1, v1);
413
414 CPPUNIT_ASSERT( all_elements(vr == y) );
415 }
416 }
417
418 template <class T>
419 void
MatrixRowVector()420 TestMathOps<T>::MatrixRowVector() {
421 vector_type r0, r1, r2;
422
423 r0 = row(m1, 0);
424 r1 = row(m1, 1);
425 r2 = row(m1, 2);
426
427 CPPUNIT_ASSERT( all_elements(r0 == m1_r0) );
428 CPPUNIT_ASSERT( all_elements(r1 == m1_r1) );
429 CPPUNIT_ASSERT( all_elements(r2 == m1_r2) );
430 }
431
432 template <class T>
433 void
MatrixColVector()434 TestMathOps<T>::MatrixColVector() {
435 vector_type c0, c1, c2;
436
437 c0 = col(m1, 0);
438 c1 = col(m1, 1);
439 c2 = col(m1, 2);
440
441 CPPUNIT_ASSERT( all_elements(c0 == m1_c0) );
442 CPPUNIT_ASSERT( all_elements(c1 == m1_c1) );
443 CPPUNIT_ASSERT( all_elements(c2 == m1_c2) );
444 }
445
446 template <class T>
447 void
MatrixDiagVector()448 TestMathOps<T>::MatrixDiagVector() {
449 vector_type vd, t;
450
451 vd = T(1), T(5), T(9);
452
453 t = diag(m1);
454
455 CPPUNIT_ASSERT( all_elements(vd == t) );
456 }
457
458 template <class T>
459 void
MatrixMatrixVector()460 TestMathOps<T>::MatrixMatrixVector() {
461 {
462 vector_type t1;
463 vector_type vr1(0), vr2(0); // clear it before use due to util::Gemv algo
464
465 // Matrix-Vector-Vector, referenz is using two ops
466 tvmet::util::Gemv(m1, v1, vr1);
467 tvmet::util::Gevvmul(vr1, vBig, vr2);
468
469 t1 = m1 * v1 * vBig;
470 CPPUNIT_ASSERT( all_elements(t1 == vr2) );
471 }
472 #if 0
473 {
474 // XXX not working due to missing operators for (XprMatrix, Vector)
475 vector_type t;
476 matrix_type vr1;
477 vector_type vr2;
478
479 // Matrix-Matrix-Vector
480 tvmet::util::Gemm(m1, mBig, vr1);
481 tvmet::util::Gemv(vr1, v1, vr2);
482
483 }
484 #endif
485 }
486
487 #endif // TVMET_TEST_MATHOPS_H
488
489 // Local Variables:
490 // mode:C++
491 // End:
492