source: ThirdParty/LinearAlgebra/src/unittests/RealSpaceMatrixUnitTest.cpp@ 9346af

Action_Thermostats ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults
Last change on this file since 9346af was 4ecb2d, checked in by Frederik Heber <heber@…>, 8 years ago

Moved LinearAlgebra sub-package into ThirdParty folder.

  • needed to adapt location of libLinearAlgebra.la in all Makefile.am's.
  • relinked m4 subfolder, relinked am_doxygen_include.am. Both point to those present in molecuilder parent folder.
  • adapted configure.ac's:
  • Property mode set to 100644
File size: 12.9 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
5 *
6 *
7 * This file is part of MoleCuilder.
8 *
9 * MoleCuilder is free software: you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation, either version 2 of the License, or
12 * (at your option) any later version.
13 *
14 * MoleCuilder is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
21 */
22
23/*
24 * RealSpaceMatrixUnitTest.cpp
25 *
26 * Created on: Jul 7, 2010
27 * Author: crueger
28 */
29
30// include config.h
31#ifdef HAVE_CONFIG_H
32#include <config.h>
33#endif
34
35#include <cppunit/CompilerOutputter.h>
36#include <cppunit/extensions/TestFactoryRegistry.h>
37#include <cppunit/ui/text/TestRunner.h>
38
39#include <cmath>
40#include <limits>
41
42#include "RealSpaceMatrixUnitTest.hpp"
43
44// include headers that implement a archive in simple text format
45#include <boost/archive/text_oarchive.hpp>
46#include <boost/archive/text_iarchive.hpp>
47
48#include "CodePatterns/Assert.hpp"
49
50#include "defs.hpp"
51#include "Exceptions.hpp"
52#include "MatrixContent.hpp"
53#include "RealSpaceMatrix.hpp"
54#include "Vector.hpp"
55
56#ifdef HAVE_TESTRUNNER
57#include "UnitTestMain.hpp"
58#endif /*HAVE_TESTRUNNER*/
59
60// Registers the fixture into the 'registry'
61CPPUNIT_TEST_SUITE_REGISTRATION( RealSpaceMatrixTest );
62
63void RealSpaceMatrixTest::setUp(){
64 // failing asserts should be thrown
65 ASSERT_DO(Assert::Throw);
66
67 zero = new RealSpaceMatrix();
68 for(int i =NDIM;i--;) {
69 for(int j =NDIM;j--;) {
70 zero->at(i,j)=0.;
71 }
72 }
73 one = new RealSpaceMatrix();
74 for(int i =NDIM;i--;){
75 one->at(i,i)=1.;
76 }
77 full=new RealSpaceMatrix();
78 for(int i=NDIM;i--;){
79 for(int j=NDIM;j--;){
80 full->at(i,j)=1.;
81 }
82 }
83 diagonal = new RealSpaceMatrix();
84 for(int i=NDIM;i--;){
85 diagonal->at(i,i)=i+1.;
86 }
87 perm1 = new RealSpaceMatrix();
88 perm1->column(0) = unitVec[0];
89 perm1->column(1) = unitVec[2];
90 perm1->column(2) = unitVec[1];
91
92
93 perm2 = new RealSpaceMatrix();
94 perm2->column(0) = unitVec[1];
95 perm2->column(1) = unitVec[0];
96 perm2->column(2) = unitVec[2];
97
98 perm3 = new RealSpaceMatrix();
99 perm3->column(0) = unitVec[1];
100 perm3->column(1) = unitVec[2];
101 perm3->column(2) = unitVec[0];
102
103 perm4 = new RealSpaceMatrix();
104 perm4->column(0) = unitVec[2];
105 perm4->column(1) = unitVec[1];
106 perm4->column(2) = unitVec[0];
107
108 perm5 = new RealSpaceMatrix();
109 perm5->column(0) = unitVec[2];
110 perm5->column(1) = unitVec[0];
111 perm5->column(2) = unitVec[1];
112
113}
114void RealSpaceMatrixTest::tearDown(){
115 delete zero;
116 delete one;
117 delete full;
118 delete diagonal;
119 delete perm1;
120 delete perm2;
121 delete perm3;
122 delete perm4;
123 delete perm5;
124}
125
126void RealSpaceMatrixTest::AccessTest(){
127 RealSpaceMatrix mat;
128 for(int i=NDIM;i--;){
129 for(int j=NDIM;j--;){
130 CPPUNIT_ASSERT_EQUAL(mat.at(i,j),0.);
131 }
132 }
133 int k=1;
134 for(int i=NDIM;i--;){
135 for(int j=NDIM;j--;){
136 mat.at(i,j)=k++;
137 }
138 }
139 k=1;
140 for(int i=NDIM;i--;){
141 for(int j=NDIM;j--;){
142 CPPUNIT_ASSERT_EQUAL(mat.at(i,j),(double)k);
143 ++k;
144 }
145 }
146}
147
148void RealSpaceMatrixTest::VectorTest(){
149 RealSpaceMatrix mat;
150 for(int i=NDIM;i--;){
151 CPPUNIT_ASSERT_EQUAL(mat.row(i),zeroVec);
152 CPPUNIT_ASSERT_EQUAL(mat.column(i),zeroVec);
153 }
154 CPPUNIT_ASSERT_EQUAL(mat.diagonal(),zeroVec);
155
156 mat.setIdentity();
157 CPPUNIT_ASSERT_EQUAL(mat.row(0),unitVec[0]);
158 CPPUNIT_ASSERT_EQUAL(mat.row(1),unitVec[1]);
159 CPPUNIT_ASSERT_EQUAL(mat.row(2),unitVec[2]);
160 CPPUNIT_ASSERT_EQUAL(mat.column(0),unitVec[0]);
161 CPPUNIT_ASSERT_EQUAL(mat.column(1),unitVec[1]);
162 CPPUNIT_ASSERT_EQUAL(mat.column(2),unitVec[2]);
163
164 Vector t1=Vector(1.,1.,1.);
165 Vector t2=Vector(2.,2.,2.);
166 Vector t3=Vector(3.,3.,3.);
167 Vector t4=Vector(1.,2.,3.);
168
169 mat.row(0)=t1;
170 mat.row(1)=t2;
171 mat.row(2)=t3;
172 CPPUNIT_ASSERT_EQUAL(mat.row(0),t1);
173 CPPUNIT_ASSERT_EQUAL(mat.row(1),t2);
174 CPPUNIT_ASSERT_EQUAL(mat.row(2),t3);
175 CPPUNIT_ASSERT_EQUAL(mat.column(0),t4);
176 CPPUNIT_ASSERT_EQUAL(mat.column(1),t4);
177 CPPUNIT_ASSERT_EQUAL(mat.column(2),t4);
178 CPPUNIT_ASSERT_EQUAL(mat.diagonal(),t4);
179 for(int i=NDIM;i--;){
180 for(int j=NDIM;j--;){
181 CPPUNIT_ASSERT_EQUAL(mat.at(i,j),i+1.);
182 }
183 }
184
185 mat.column(0)=t1;
186 mat.column(1)=t2;
187 mat.column(2)=t3;
188 CPPUNIT_ASSERT_EQUAL(mat.column(0),t1);
189 CPPUNIT_ASSERT_EQUAL(mat.column(1),t2);
190 CPPUNIT_ASSERT_EQUAL(mat.column(2),t3);
191 CPPUNIT_ASSERT_EQUAL(mat.row(0),t4);
192 CPPUNIT_ASSERT_EQUAL(mat.row(1),t4);
193 CPPUNIT_ASSERT_EQUAL(mat.row(2),t4);
194 CPPUNIT_ASSERT_EQUAL(mat.diagonal(),t4);
195 for(int i=NDIM;i--;){
196 for(int j=NDIM;j--;){
197 CPPUNIT_ASSERT_EQUAL(mat.at(i,j),j+1.);
198 }
199 }
200}
201
202void RealSpaceMatrixTest::TransposeTest(){
203 RealSpaceMatrix res;
204
205 // transpose of unit is unit
206 res.setIdentity();
207 res.transpose();
208 CPPUNIT_ASSERT_EQUAL(res,*one);
209
210 // transpose of transpose is same matrix
211 res.setZero();
212 res.set(2,2, 1.);
213 CPPUNIT_ASSERT_EQUAL(res.transpose().transpose(),res);
214}
215
216void RealSpaceMatrixTest::OperationTest(){
217 RealSpaceMatrix res;
218
219 res =(*zero) *(*zero);
220 //std::cout << *zero << " times " << *zero << " is " << res << std::endl;
221 CPPUNIT_ASSERT_EQUAL(res,*zero);
222 res =(*zero) *(*one);
223 CPPUNIT_ASSERT_EQUAL(res,*zero);
224 res =(*zero) *(*full);
225 CPPUNIT_ASSERT_EQUAL(res,*zero);
226 res =(*zero) *(*diagonal);
227 CPPUNIT_ASSERT_EQUAL(res,*zero);
228 res =(*zero) *(*perm1);
229 CPPUNIT_ASSERT_EQUAL(res,*zero);
230 res =(*zero) *(*perm2);
231 CPPUNIT_ASSERT_EQUAL(res,*zero);
232 res =(*zero) *(*perm3);
233 CPPUNIT_ASSERT_EQUAL(res,*zero);
234 res =(*zero) *(*perm4);
235 CPPUNIT_ASSERT_EQUAL(res,*zero);
236 res =(*zero) *(*perm5);
237 CPPUNIT_ASSERT_EQUAL(res,*zero);
238
239 res =(*one)*(*one);
240 CPPUNIT_ASSERT_EQUAL(res,*one);
241 res =(*one)*(*full);
242 CPPUNIT_ASSERT_EQUAL(res,*full);
243 res =(*one)*(*diagonal);
244 CPPUNIT_ASSERT_EQUAL(res,*diagonal);
245 res =(*one)*(*perm1);
246 CPPUNIT_ASSERT_EQUAL(res,*perm1);
247 res =(*one)*(*perm2);
248 CPPUNIT_ASSERT_EQUAL(res,*perm2);
249 res =(*one)*(*perm3);
250 CPPUNIT_ASSERT_EQUAL(res,*perm3);
251 res =(*one)*(*perm4);
252 CPPUNIT_ASSERT_EQUAL(res,*perm4);
253 res =(*one)*(*perm5);
254 CPPUNIT_ASSERT_EQUAL(res,*perm5);
255
256 res = (*full)*(*perm1);
257 CPPUNIT_ASSERT_EQUAL(res,*full);
258 res = (*full)*(*perm2);
259 CPPUNIT_ASSERT_EQUAL(res,*full);
260 res = (*full)*(*perm3);
261 CPPUNIT_ASSERT_EQUAL(res,*full);
262 res = (*full)*(*perm4);
263 CPPUNIT_ASSERT_EQUAL(res,*full);
264 res = (*full)*(*perm5);
265 CPPUNIT_ASSERT_EQUAL(res,*full);
266
267 res = (*diagonal)*(*perm1);
268 CPPUNIT_ASSERT_EQUAL(res.column(0),unitVec[0]);
269 CPPUNIT_ASSERT_EQUAL(res.column(1),3*unitVec[2]);
270 CPPUNIT_ASSERT_EQUAL(res.column(2),2*unitVec[1]);
271 res = (*diagonal)*(*perm2);
272 CPPUNIT_ASSERT_EQUAL(res.column(0),2*unitVec[1]);
273 CPPUNIT_ASSERT_EQUAL(res.column(1),unitVec[0]);
274 CPPUNIT_ASSERT_EQUAL(res.column(2),3*unitVec[2]);
275 res = (*diagonal)*(*perm3);
276 CPPUNIT_ASSERT_EQUAL(res.column(0),2*unitVec[1]);
277 CPPUNIT_ASSERT_EQUAL(res.column(1),3*unitVec[2]);
278 CPPUNIT_ASSERT_EQUAL(res.column(2),unitVec[0]);
279 res = (*diagonal)*(*perm4);
280 CPPUNIT_ASSERT_EQUAL(res.column(0),3*unitVec[2]);
281 CPPUNIT_ASSERT_EQUAL(res.column(1),2*unitVec[1]);
282 CPPUNIT_ASSERT_EQUAL(res.column(2),unitVec[0]);
283 res = (*diagonal)*(*perm5);
284 CPPUNIT_ASSERT_EQUAL(res.column(0),3*unitVec[2]);
285 CPPUNIT_ASSERT_EQUAL(res.column(1),unitVec[0]);
286 CPPUNIT_ASSERT_EQUAL(res.column(2),2*unitVec[1]);
287}
288
289void RealSpaceMatrixTest::RotationTest(){
290 RealSpaceMatrix res;
291 RealSpaceMatrix inverse;
292 Vector fixture;
293
294 // zero rotation angles yields unity matrix
295 res.setRotation(0,0,0);
296 CPPUNIT_ASSERT_EQUAL(*one, res);
297
298 // arbitrary rotation matrix has det = 1
299 res.setRotation(M_PI/3.,1.,M_PI/7.);
300 CPPUNIT_ASSERT(fabs(res.determinant() - 1.) <= LINALG_MYEPSILON());
301
302 // inverse is rotation matrix with negative angles
303 res.setRotation(M_PI/3.,0.,0.);
304 inverse.setRotation(-M_PI/3.,0.,0.);
305 CPPUNIT_ASSERT_EQUAL(*one, res * inverse);
306
307 // ... or transposed
308 res.setRotation(M_PI/3.,0.,0.);
309 CPPUNIT_ASSERT_EQUAL(inverse, res.transposed());
310
311 // check arbitrary axis
312 res.setRotation(unitVec[0], M_PI/3.);
313 inverse.setRotation(M_PI/3.,0.,0.);
314 CPPUNIT_ASSERT_EQUAL( res, inverse);
315 res.setRotation(unitVec[1], M_PI/3.);
316 inverse.setRotation(0., M_PI/3.,0.);
317 CPPUNIT_ASSERT_EQUAL( res, inverse);
318 res.setRotation(unitVec[2], M_PI/3.);
319 inverse.setRotation(0.,0.,M_PI/3.);
320 CPPUNIT_ASSERT_EQUAL( res, inverse);
321
322 // check axis not normalized throws
323#ifndef NDEBUG
324 CPPUNIT_ASSERT_THROW( res.setRotation(unitVec[0]+unitVec[1], M_PI/3.),Assert::AssertionFailure);
325#endif
326
327 // check arbitrary axis and det=1
328 fixture = (unitVec[0]+unitVec[1]).getNormalized();
329 res.setRotation(fixture, M_PI/3.);
330 CPPUNIT_ASSERT(fabs(res.determinant() - 1.) <= LINALG_MYEPSILON());
331 fixture = (unitVec[1]+unitVec[2]).getNormalized();
332 res.setRotation(fixture, M_PI/3.);
333 CPPUNIT_ASSERT(fabs(res.determinant() - 1.) <= LINALG_MYEPSILON());
334 fixture = (unitVec[1]+unitVec[2]).getNormalized();
335 res.setRotation(fixture, M_PI/3.);
336 CPPUNIT_ASSERT(fabs(res.determinant() - 1.) <= LINALG_MYEPSILON());
337}
338
339void RealSpaceMatrixTest::InvertTest(){
340 CPPUNIT_ASSERT_THROW(zero->invert(),NotInvertibleException);
341 CPPUNIT_ASSERT_THROW(full->invert(),NotInvertibleException);
342
343 RealSpaceMatrix res;
344 res = (*one)*one->invert();
345 CPPUNIT_ASSERT_EQUAL(res,*one);
346 res = (*diagonal)*diagonal->invert();
347 CPPUNIT_ASSERT_EQUAL(res,*one);
348 res = (*perm1)*perm1->invert();
349 CPPUNIT_ASSERT_EQUAL(res,*one);
350 res = (*perm2)*perm2->invert();
351 CPPUNIT_ASSERT_EQUAL(res,*one);
352 res = (*perm3)*perm3->invert();
353 CPPUNIT_ASSERT_EQUAL(res,*one);
354 res = (*perm4)*perm4->invert();
355 CPPUNIT_ASSERT_EQUAL(res,*one);
356 res = (*perm5)*perm5->invert();
357 CPPUNIT_ASSERT_EQUAL(res,*one);
358}
359
360
361void RealSpaceMatrixTest::DeterminantTest(){
362 CPPUNIT_ASSERT_EQUAL(zero->determinant(),0.);
363 CPPUNIT_ASSERT_EQUAL(one->determinant(),1.);
364 CPPUNIT_ASSERT_EQUAL(diagonal->determinant(),6.);
365 CPPUNIT_ASSERT_EQUAL(full->determinant(),0.);
366 CPPUNIT_ASSERT_EQUAL(perm1->determinant(),-1.);
367 CPPUNIT_ASSERT_EQUAL(perm2->determinant(),-1.);
368 CPPUNIT_ASSERT_EQUAL(perm3->determinant(),1.);
369 CPPUNIT_ASSERT_EQUAL(perm4->determinant(),-1.);
370 CPPUNIT_ASSERT_EQUAL(perm5->determinant(),1.);
371}
372
373void RealSpaceMatrixTest::VecMultTest(){
374 CPPUNIT_ASSERT_EQUAL((*zero)*unitVec[0],zeroVec);
375 CPPUNIT_ASSERT_EQUAL((*zero)*unitVec[1],zeroVec);
376 CPPUNIT_ASSERT_EQUAL((*zero)*unitVec[2],zeroVec);
377 CPPUNIT_ASSERT_EQUAL((*zero)*zeroVec,zeroVec);
378
379 CPPUNIT_ASSERT_EQUAL((*one)*unitVec[0],unitVec[0]);
380 CPPUNIT_ASSERT_EQUAL((*one)*unitVec[1],unitVec[1]);
381 CPPUNIT_ASSERT_EQUAL((*one)*unitVec[2],unitVec[2]);
382 CPPUNIT_ASSERT_EQUAL((*one)*zeroVec,zeroVec);
383
384 CPPUNIT_ASSERT_EQUAL((*diagonal)*unitVec[0],unitVec[0]);
385 CPPUNIT_ASSERT_EQUAL((*diagonal)*unitVec[1],2*unitVec[1]);
386 CPPUNIT_ASSERT_EQUAL((*diagonal)*unitVec[2],3*unitVec[2]);
387 CPPUNIT_ASSERT_EQUAL((*diagonal)*zeroVec,zeroVec);
388
389 CPPUNIT_ASSERT_EQUAL((*perm1)*unitVec[0],unitVec[0]);
390 CPPUNIT_ASSERT_EQUAL((*perm1)*unitVec[1],unitVec[2]);
391 CPPUNIT_ASSERT_EQUAL((*perm1)*unitVec[2],unitVec[1]);
392 CPPUNIT_ASSERT_EQUAL((*perm1)*zeroVec,zeroVec);
393
394 CPPUNIT_ASSERT_EQUAL((*perm2)*unitVec[0],unitVec[1]);
395 CPPUNIT_ASSERT_EQUAL((*perm2)*unitVec[1],unitVec[0]);
396 CPPUNIT_ASSERT_EQUAL((*perm2)*unitVec[2],unitVec[2]);
397 CPPUNIT_ASSERT_EQUAL((*perm2)*zeroVec,zeroVec);
398
399 CPPUNIT_ASSERT_EQUAL((*perm3)*unitVec[0],unitVec[1]);
400 CPPUNIT_ASSERT_EQUAL((*perm3)*unitVec[1],unitVec[2]);
401 CPPUNIT_ASSERT_EQUAL((*perm3)*unitVec[2],unitVec[0]);
402 CPPUNIT_ASSERT_EQUAL((*perm3)*zeroVec,zeroVec);
403
404 CPPUNIT_ASSERT_EQUAL((*perm4)*unitVec[0],unitVec[2]);
405 CPPUNIT_ASSERT_EQUAL((*perm4)*unitVec[1],unitVec[1]);
406 CPPUNIT_ASSERT_EQUAL((*perm4)*unitVec[2],unitVec[0]);
407 CPPUNIT_ASSERT_EQUAL((*perm4)*zeroVec,zeroVec);
408
409 CPPUNIT_ASSERT_EQUAL((*perm5)*unitVec[0],unitVec[2]);
410 CPPUNIT_ASSERT_EQUAL((*perm5)*unitVec[1],unitVec[0]);
411 CPPUNIT_ASSERT_EQUAL((*perm5)*unitVec[2],unitVec[1]);
412 CPPUNIT_ASSERT_EQUAL((*perm5)*zeroVec,zeroVec);
413
414 Vector t = Vector(1.,2.,3.);
415 CPPUNIT_ASSERT_EQUAL((*perm1)*t,Vector(1,3,2));
416 CPPUNIT_ASSERT_EQUAL((*perm2)*t,Vector(2,1,3));
417 CPPUNIT_ASSERT_EQUAL((*perm3)*t,Vector(3,1,2));
418 CPPUNIT_ASSERT_EQUAL((*perm4)*t,Vector(3,2,1));
419 CPPUNIT_ASSERT_EQUAL((*perm5)*t,Vector(2,3,1));
420}
421
422void RealSpaceMatrixTest::SerializationTest()
423{
424 // write element to stream
425 std::stringstream stream;
426 boost::archive::text_oarchive oa(stream);
427 oa << diagonal;
428
429 //std::cout << "Contents of archive is " << stream.str() << std::endl;
430
431 // create and open an archive for input
432 boost::archive::text_iarchive ia(stream);
433 // read class state from archive
434 RealSpaceMatrix *newm;
435
436 ia >> newm;
437
438 CPPUNIT_ASSERT (*diagonal == *newm);
439}
Note: See TracBrowser for help on using the repository browser.