OpenWalnut 1.2.5
|
00001 //--------------------------------------------------------------------------- 00002 // 00003 // Project: OpenWalnut ( http://www.openwalnut.org ) 00004 // 00005 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS 00006 // For more information see http://www.openwalnut.org/copying 00007 // 00008 // This file is part of OpenWalnut. 00009 // 00010 // OpenWalnut is free software: you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as published by 00012 // the Free Software Foundation, either version 3 of the License, or 00013 // (at your option) any later version. 00014 // 00015 // OpenWalnut is distributed in the hope that it will be useful, 00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00018 // GNU Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public License 00021 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>. 00022 // 00023 //--------------------------------------------------------------------------- 00024 00025 #ifndef WTENSORFUNCTIONS_TEST_H 00026 #define WTENSORFUNCTIONS_TEST_H 00027 00028 #include <string> 00029 #include <vector> 00030 #include <algorithm> 00031 00032 #include <cxxtest/TestSuite.h> 00033 00034 #include "../../WAssert.h" 00035 #include "../../WLimits.h" 00036 #include "../../WStringUtils.h" 00037 #include "../WTensorFunctions.h" 00038 00039 /** 00040 * Test class for some tensor functions. 00041 */ 00042 class WTensorFunctionsTest : public CxxTest::TestSuite 00043 { 00044 public: 00045 /** 00046 * The eigenvalue of the symmetrical matrix: 00047 * 0.000179516, 2.09569e-05, 2.76557e-06, 0.000170189, -5.52619e-07, 0.00015239 00048 * (0.000196397;0.000155074;0.000150625) 00049 */ 00050 void testSpecialSymMatrixEigenvalueTestCaseNumericalStability( void ) 00051 { 00052 WTensorSym< 2, 3 > t; 00053 t( 0, 0 ) = 0.000179516; 00054 t( 0, 1 ) = 2.09569e-05; 00055 t( 0, 2 ) = 2.76557e-06; 00056 t( 1, 1 ) = 0.000170189; 00057 t( 1, 2 ) = -5.52619e-07; 00058 t( 2, 2 ) = 0.00015239; 00059 RealEigenSystem sys; 00060 jacobiEigenvector3D( t, &sys ); 00061 00062 TS_ASSERT_DELTA( sys[0].first, 1.5062467240725114e-04, 1e-9 ); 00063 TS_ASSERT_DELTA( sys[1].first, 1.5507354000104679e-04, 1e-9 ); 00064 TS_ASSERT_DELTA( sys[2].first, 1.9639678759170208e-04, 1e-9 ); 00065 } 00066 00067 /** 00068 * Test the jacobi eigenvector calculation. 00069 */ 00070 void testJacobiEigenvectors() 00071 { 00072 // the test matrix 00073 WTensorSym< 2, 3 > t; 00074 00075 // simple diagonal matrices should be decomposed correctly 00076 // 1 2 3 00077 t( 0, 0 ) = 1.0; 00078 t( 1, 1 ) = 2.0; 00079 t( 2, 2 ) = 3.0; 00080 00081 RealEigenSystem sys; 00082 jacobiEigenvector3D( t, &sys ); 00083 00084 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 ); 00085 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 ); 00086 TS_ASSERT_DELTA( sys[2].first, 3.0, 1e-6 ); 00087 // eigenvectors should be perpendicular 00088 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00089 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00090 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00091 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00092 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00093 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00094 00095 // 1 2 -3 00096 t( 0, 0 ) = 1.0; 00097 t( 1, 1 ) = 2.0; 00098 t( 2, 2 ) = -3.0; 00099 00100 jacobiEigenvector3D( t, &sys ); 00101 00102 TS_ASSERT_DELTA( sys[0].first, -3.0, 1e-6 ); 00103 TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 ); 00104 TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 ); 00105 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00106 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00107 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00108 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00109 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00110 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00111 00112 // 1 2 2 00113 t( 0, 0 ) = 1.0; 00114 t( 1, 1 ) = 2.0; 00115 t( 2, 2 ) = 2.0; 00116 00117 jacobiEigenvector3D( t, &sys ); 00118 00119 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 ); 00120 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 ); 00121 TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 ); 00122 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00123 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00124 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00125 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00126 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00127 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00128 00129 // -1 -1 -1 00130 t( 0, 0 ) = -1.0; 00131 t( 1, 1 ) = -1.0; 00132 t( 2, 2 ) = -1.0; 00133 00134 jacobiEigenvector3D( t, &sys ); 00135 00136 TS_ASSERT_DELTA( sys[0].first, -1.0, 1e-6 ); 00137 TS_ASSERT_DELTA( sys[1].first, -1.0, 1e-6 ); 00138 TS_ASSERT_DELTA( sys[2].first, -1.0, 1e-6 ); 00139 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00140 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00141 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00142 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00143 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00144 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00145 00146 // 1 0 1 00147 t( 0, 0 ) = 1.0; 00148 t( 1, 1 ) = 0.0; 00149 t( 2, 2 ) = 1.0; 00150 00151 jacobiEigenvector3D( t, &sys ); 00152 00153 TS_ASSERT_DELTA( sys[0].first, 0.0, 1e-6 ); 00154 TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 ); 00155 TS_ASSERT_DELTA( sys[2].first, 1.0, 1e-6 ); 00156 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00157 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00158 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00159 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00160 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00161 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00162 00163 // 0 0 0 00164 t( 0, 0 ) = 0.0; 00165 t( 1, 1 ) = 0.0; 00166 t( 2, 2 ) = 0.0; 00167 00168 jacobiEigenvector3D( t, &sys ); 00169 00170 TS_ASSERT_DELTA( sys[0].first, 0.0, 1e-6 ); 00171 TS_ASSERT_DELTA( sys[1].first, 0.0, 1e-6 ); 00172 TS_ASSERT_DELTA( sys[2].first, 0.0, 1e-6 ); 00173 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00174 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00175 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00176 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00177 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00178 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00179 00180 // similar eigenvalues 00181 // 2.000001 0.0 1.999998 00182 t( 0, 0 ) = 2.000001; 00183 t( 1, 1 ) = 0.0; 00184 t( 2, 2 ) = 1.999998; 00185 00186 jacobiEigenvector3D( t, &sys ); 00187 00188 TS_ASSERT_DELTA( sys[0].first, 0.0, 1e-6 ); 00189 TS_ASSERT_DELTA( sys[1].first, 1.999998, 1e-6 ); 00190 TS_ASSERT_DELTA( sys[2].first, 2.000001, 1e-6 ); 00191 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00192 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00193 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00194 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00195 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00196 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00197 00198 // very large eigenvalues 00199 // 3.824572321236e1000 1 2 00200 t( 0, 0 ) = 3.824572321236e30; 00201 t( 1, 1 ) = 1.0; 00202 t( 2, 2 ) = 2.0; 00203 00204 jacobiEigenvector3D( t, &sys ); 00205 00206 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 ); 00207 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 ); 00208 TS_ASSERT_DELTA( sys[2].first, 3.824572321236e30, 1e-6 ); 00209 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00210 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00211 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00212 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00213 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00214 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00215 00216 // very small eigenvalues 00217 // 3.824572321236e-1000 1 2 00218 t( 0, 0 ) = 3.824572321236e-30; 00219 t( 1, 1 ) = 1.0; 00220 t( 2, 2 ) = 2.0; 00221 00222 jacobiEigenvector3D( t, &sys ); 00223 00224 TS_ASSERT_DELTA( sys[0].first, 3.824572321236e-30, 1e-6 ); 00225 TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 ); 00226 TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 ); 00227 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00228 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00229 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00230 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00231 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00232 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00233 00234 // some more sophisticated tests 00235 // (using similarity transformations on diagonal matrices to create test cases) 00236 t( 0, 0 ) = 1; 00237 t( 1, 1 ) = 2; 00238 t( 2, 2 ) = 3; 00239 00240 // note that the jacobi eigenvector functions does not sort the eigenvalues and 00241 // eigenvectors that were found 00242 t = similarity_rotate_givens( t, 0, 2, 2.78 ); 00243 00244 jacobiEigenvector3D( t, &sys ); 00245 00246 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 ); 00247 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 ); 00248 TS_ASSERT_DELTA( sys[2].first, 3.0, 1e-6 ); 00249 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00250 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00251 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00252 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00253 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00254 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00255 compare_results( t, sys ); 00256 00257 t = WTensorSym< 2, 3 >(); 00258 t( 0, 0 ) = 2; 00259 t( 1, 1 ) = 1; 00260 t( 2, 2 ) = 3; 00261 00262 t = similarity_rotate_givens( t, 0, 2, 2.79 ); 00263 t = similarity_rotate_givens( t, 1, 2, -3.44 ); 00264 t = similarity_rotate_givens( t, 1, 0, -0.46 ); 00265 t = similarity_rotate_givens( t, 2, 1, 5.98 ); 00266 00267 jacobiEigenvector3D( t, &sys ); 00268 00269 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 ); 00270 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 ); 00271 TS_ASSERT_DELTA( sys[2].first, 3.0, 1e-6 ); 00272 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00273 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00274 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00275 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00276 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00277 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00278 compare_results( t, sys ); 00279 00280 t = WTensorSym< 2, 3 >(); 00281 t( 0, 0 ) = 2; 00282 t( 1, 1 ) = 2; 00283 t( 2, 2 ) = 2; 00284 00285 t = similarity_rotate_givens( t, 1, 2, -3.44 ); 00286 t = similarity_rotate_givens( t, 2, 1, 5.98 ); 00287 t = similarity_rotate_givens( t, 0, 2, 2.79 ); 00288 t = similarity_rotate_givens( t, 1, 0, -0.46 ); 00289 00290 jacobiEigenvector3D( t, &sys ); 00291 00292 TS_ASSERT_DELTA( sys[0].first, 2.0, 1e-6 ); 00293 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 ); 00294 TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 ); 00295 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00296 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00297 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00298 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00299 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00300 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00301 compare_results( t, sys ); 00302 } 00303 00304 /** 00305 * Test the cardano eigenvalue calculation. 00306 */ 00307 void testCardanoEigenvalues() 00308 { 00309 // the test matrix 00310 WTensorSym< 2, 3 > t; 00311 // variables for the output values 00312 std::vector< double > d( 3 ); 00313 00314 // simple diagonal matrices should be decomposed correctly 00315 00316 // 1 2 3 00317 t( 0, 0 ) = 1.0; 00318 t( 1, 1 ) = 2.0; 00319 t( 2, 2 ) = 3.0; 00320 00321 d = getEigenvaluesCardano( t ); 00322 00323 TS_ASSERT_DELTA( d[ 0 ], 3.0, 1e-6 ); 00324 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 ); 00325 TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 ); 00326 00327 // 1 2 -3 00328 t( 0, 0 ) = 1.0; 00329 t( 1, 1 ) = 2.0; 00330 t( 2, 2 ) = -3.0; 00331 00332 d = getEigenvaluesCardano( t ); 00333 00334 TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 ); 00335 TS_ASSERT_DELTA( d[ 1 ], 1.0, 1e-6 ); 00336 TS_ASSERT_DELTA( d[ 2 ], -3.0, 1e-6 ); 00337 00338 // 1 2 2 00339 t( 0, 0 ) = 1.0; 00340 t( 1, 1 ) = 2.0; 00341 t( 2, 2 ) = 2.0; 00342 00343 d = getEigenvaluesCardano( t ); 00344 00345 TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 ); 00346 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 ); 00347 TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 ); 00348 00349 // -1 -1 -1 00350 t( 0, 0 ) = -1.0; 00351 t( 1, 1 ) = -1.0; 00352 t( 2, 2 ) = -1.0; 00353 00354 d = getEigenvaluesCardano( t ); 00355 00356 TS_ASSERT_DELTA( d[ 0 ], -1.0, 1e-6 ); 00357 TS_ASSERT_DELTA( d[ 1 ], -1.0, 1e-6 ); 00358 TS_ASSERT_DELTA( d[ 2 ], -1.0, 1e-6 ); 00359 00360 // 1 0 1 00361 t( 0, 0 ) = 1.0; 00362 t( 1, 1 ) = 0.0; 00363 t( 2, 2 ) = 1.0; 00364 00365 d = getEigenvaluesCardano( t ); 00366 00367 TS_ASSERT_DELTA( d[ 0 ], 1.0, 1e-6 ); 00368 TS_ASSERT_DELTA( d[ 1 ], 1.0, 1e-6 ); 00369 TS_ASSERT_DELTA( d[ 2 ], 0.0, 1e-6 ); 00370 00371 // 0 0 0 00372 t( 0, 0 ) = 0.0; 00373 t( 1, 1 ) = 0.0; 00374 t( 2, 2 ) = 0.0; 00375 00376 d = getEigenvaluesCardano( t ); 00377 00378 TS_ASSERT_DELTA( d[ 0 ], 0.0, 1e-6 ); 00379 TS_ASSERT_DELTA( d[ 1 ], 0.0, 1e-6 ); 00380 TS_ASSERT_DELTA( d[ 2 ], 0.0, 1e-6 ); 00381 00382 // similar eigenvalues 00383 // 2.000001 0.0 1.999998 00384 t( 0, 0 ) = 2.000001; 00385 t( 1, 1 ) = 0.0; 00386 t( 2, 2 ) = 1.999998; 00387 00388 d = getEigenvaluesCardano( t ); 00389 00390 TS_ASSERT_DELTA( d[ 0 ], 2.000001, 1e-6 ); 00391 TS_ASSERT_DELTA( d[ 1 ], 1.999998, 1e-6 ); 00392 TS_ASSERT_DELTA( d[ 2 ], 0.0, 1e-6 ); 00393 00394 // very large eigenvalues 00395 // 3.824572321236e10 1 2 00396 t( 0, 0 ) = 3.824572321236e10; 00397 t( 1, 1 ) = 1.0; 00398 t( 2, 2 ) = 2.0; 00399 00400 d = getEigenvaluesCardano( t ); 00401 00402 TS_ASSERT_DELTA( d[ 0 ], 3.824572321236e10, 1e-6 ); 00403 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 ); 00404 TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 ); 00405 00406 // very small eigenvalues 00407 // 3.824572321236e-10 1 2 00408 t( 0, 0 ) = 3.824572321236e-30; 00409 t( 1, 1 ) = 1.0; 00410 t( 2, 2 ) = 2.0; 00411 00412 d = getEigenvaluesCardano( t ); 00413 00414 TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 ); 00415 TS_ASSERT_DELTA( d[ 1 ], 1.0, 1e-6 ); 00416 TS_ASSERT_DELTA( d[ 2 ], 3.824572321236e-30, 1e-6 ); 00417 00418 // some more sophisticated tests 00419 // (using similarity transformations on diagonal matrices to create test cases) 00420 t( 0, 0 ) = 1; 00421 t( 1, 1 ) = 2; 00422 t( 2, 2 ) = 3; 00423 00424 // note that the jacobi eigenvector functions does not sort the eigenvalues and 00425 // eigenvectors that were found 00426 t = similarity_rotate_givens( t, 0, 2, 2.78 ); 00427 00428 d = getEigenvaluesCardano( t ); 00429 00430 TS_ASSERT_DELTA( d[ 0 ], 3.0, 1e-6 ); 00431 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 ); 00432 TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 ); 00433 00434 // the eigenvalues calculated for this test are a bit off 00435 //t = WTensorSym< 2, 3 >(); 00436 //t( 0, 0 ) = 2; 00437 //t( 1, 1 ) = 1; 00438 //t( 2, 2 ) = 3; 00439 00440 //t = similarity_rotate_givens( t, 0, 2, 2.79 ); 00441 //t = similarity_rotate_givens( t, 1, 2, -3.44 ); 00442 //t = similarity_rotate_givens( t, 1, 0, -0.46 ); 00443 //t = similarity_rotate_givens( t, 2, 1, 5.98 ); 00444 00445 //d = getEigenvaluesCardano( t ); 00446 00447 //TS_ASSERT_DELTA( d[ 0 ], 3.0, 1e-6 ); 00448 //TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 ); 00449 //TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 ); 00450 00451 t = WTensorSym< 2, 3 >(); 00452 t( 0, 0 ) = 2; 00453 t( 1, 1 ) = 2; 00454 t( 2, 2 ) = 2; 00455 00456 t = similarity_rotate_givens( t, 1, 2, -3.44 ); 00457 t = similarity_rotate_givens( t, 2, 1, 5.98 ); 00458 t = similarity_rotate_givens( t, 0, 2, 2.79 ); 00459 t = similarity_rotate_givens( t, 1, 0, -0.46 ); 00460 00461 d = getEigenvaluesCardano( t ); 00462 00463 TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 ); 00464 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 ); 00465 TS_ASSERT_DELTA( d[ 2 ], 2.0, 1e-6 ); 00466 } 00467 00468 private: 00469 /** 00470 * A helper function performing a similarity transform using a givens rotation. 00471 * 00472 * \param m The symmetric tensor to transform. 00473 * \param i A row index. 00474 * \param j A column index. 00475 * \param angle The rotation angle (in radians). 00476 * 00477 * \note i must not have the same values as j 00478 * 00479 * \return The new tensor 00480 */ 00481 template< std::size_t dim, typename Data_T > 00482 WTensorSym< 2, dim, Data_T > similarity_rotate_givens( WTensorSym< 2, dim, Data_T > const& m, 00483 std::size_t i, 00484 std::size_t j, 00485 double angle ) 00486 { 00487 WAssert( i != j, "" ); 00488 00489 double s = sin( angle ); 00490 double c = cos( angle ); 00491 WTensor< 2, dim, Data_T > r; 00492 for( std::size_t k = 0; k < dim; ++k ) 00493 { 00494 if( k == i || k == j ) 00495 { 00496 r( k, k ) = c; 00497 } 00498 else 00499 { 00500 r( k, k ) = 1.0; 00501 } 00502 } 00503 r( i, j ) = s; 00504 r( j, i ) = -s; 00505 WTensor< 2, dim, Data_T > t( r ); 00506 std::swap( t( i, j ), t( j, i ) ); 00507 00508 r = t * m * r; 00509 WTensorSym< 2, dim, Data_T > res; 00510 res( 0, 0 ) = r( 0, 0 ); 00511 res( 1, 0 ) = r( 1, 0 ); 00512 res( 2, 0 ) = r( 2, 0 ); 00513 res( 1, 1 ) = r( 1, 1 ); 00514 res( 2, 1 ) = r( 2, 1 ); 00515 res( 2, 2 ) = r( 2, 2 ); 00516 00517 return res; 00518 } 00519 00520 /** 00521 * Test if the given vectors are eigenvectors to the given eigenvalues of a 00522 * symmetric matrix. 00523 * 00524 * \param m A symmetric matrix. 00525 * \param sys The eigen system ( eigenvalues and eigenvectors ) 00526 */ 00527 template< std::size_t dim, typename Data_T > 00528 void compare_results( WTensorSym< 2, dim, Data_T > const& m, RealEigenSystem const& sys ) 00529 { 00530 WTensor< 2, dim, Data_T > t, r; 00531 for( std::size_t i = 0; i < dim; ++i ) 00532 { 00533 for( std::size_t j = 0; j < dim; ++j ) 00534 { 00535 t( j, i ) = sys[i].second[j]; 00536 r( j, i ) = sys[i].first * sys[i].second[j]; 00537 } 00538 } 00539 00540 t = m * t; 00541 for( std::size_t i = 0; i < dim; ++i ) 00542 { 00543 for( std::size_t j = 0; j < dim; ++j ) 00544 { 00545 TS_ASSERT_DELTA( t( i, j ), r( i, j ), 1e-15 ); 00546 } 00547 } 00548 } 00549 }; 00550 00551 /** 00552 * Test class for all tensor operators. 00553 */ 00554 class WTensorOperatorsTest : public CxxTest::TestSuite 00555 { 00556 public: 00557 /** 00558 * Test order 2 tensor multiplication. 00559 */ 00560 void testMultiplyTensorsOperator() 00561 { 00562 // some special cases, WTensor * WTensor 00563 TS_ASSERT_EQUALS( zero * zero, zero ); 00564 TS_ASSERT_EQUALS( zero * one, zero ); 00565 TS_ASSERT_EQUALS( one * zero, zero ); 00566 TS_ASSERT_EQUALS( one * one, one ); 00567 TS_ASSERT_EQUALS( one * rdm1, rdm1 ); 00568 TS_ASSERT_EQUALS( one * rdm2, rdm2 ); 00569 TS_ASSERT_EQUALS( rdm1 * one, rdm1 ); 00570 TS_ASSERT_EQUALS( rdm2 * one, rdm2 ); 00571 00572 // a normal case 00573 TS_ASSERT_EQUALS( rdm1 * rdm2, res1 ); 00574 00575 // some special cases, WTensorSym * WTensorSym 00576 TS_ASSERT_EQUALS( szero * szero, zero ); 00577 TS_ASSERT_EQUALS( szero * sone, zero ); 00578 TS_ASSERT_EQUALS( sone * szero, zero ); 00579 TS_ASSERT_EQUALS( sone * sone, one ); 00580 TS_ASSERT_EQUALS( sone * srdm1, res2 ); 00581 TS_ASSERT_EQUALS( srdm1 * sone, res2 ); 00582 TS_ASSERT_EQUALS( srdm1 * srdm2, res3 ); 00583 00584 // WTensor * WTensorSym 00585 TS_ASSERT_EQUALS( zero * szero, zero ); 00586 TS_ASSERT_EQUALS( one * sone, one ); 00587 TS_ASSERT_EQUALS( zero * sone, zero ); 00588 TS_ASSERT_EQUALS( one * szero, zero ); 00589 00590 // WTensorSym * WTensor 00591 TS_ASSERT_EQUALS( szero * zero, zero ); 00592 TS_ASSERT_EQUALS( sone * one, one ); 00593 TS_ASSERT_EQUALS( szero * one, zero ); 00594 TS_ASSERT_EQUALS( sone * zero, zero ); 00595 TS_ASSERT_EQUALS( srdm1 * rdm1, res4 ); 00596 } 00597 00598 /** 00599 * The optimizations for symmetric tensors should not corrupt the result. 00600 */ 00601 void testEvaluateSphericalFunction() 00602 { 00603 WTensorSym< 4, 3, double > t; 00604 // the tensor 00605 t( 0, 0, 0, 0 ) = 2.5476; 00606 t( 1, 1, 1, 1 ) = 3.5476; 00607 t( 2, 2, 2, 2 ) = 4.5476; 00608 t( 0, 0, 0, 1 ) = 5.5476; 00609 t( 0, 0, 0, 2 ) = 6.5476; 00610 t( 1, 1, 1, 0 ) = 7.5476; 00611 t( 1, 1, 1, 2 ) = 8.5476; 00612 t( 2, 2, 2, 0 ) = 9.5476; 00613 t( 2, 2, 2, 1 ) = 10.5476; 00614 t( 0, 0, 1, 2 ) = 11.5476; 00615 t( 1, 1, 0, 2 ) = 12.5476; 00616 t( 2, 2, 0, 1 ) = 13.5476; 00617 t( 0, 0, 1, 1 ) = 14.5476; 00618 t( 0, 0, 2, 2 ) = 15.5476; 00619 t( 1, 1, 2, 2 ) = 16.5476; 00620 00621 // the gradients 00622 std::vector< WVector3d > gradients; 00623 gradients.push_back( WVector3d( 1.0, 0.0, 0.0 ) ); 00624 gradients.push_back( WVector3d( 0.0, 1.0, 0.0 ) ); 00625 gradients.push_back( normalize( WVector3d( 1.0, 1.0, 0.0 ) ) ); 00626 gradients.push_back( normalize( WVector3d( 0.3, 0.4, 0.5 ) ) ); 00627 gradients.push_back( normalize( WVector3d( -7.0, 3.0, -1.0 ) ) ); 00628 00629 for( int k = 0; k < 5; ++k ) 00630 { 00631 double res = calcTens( t, gradients[ k ] ); 00632 TS_ASSERT_DELTA( res, evaluateSphericalFunction( t, gradients[ k ] ), 0.001 ); 00633 } 00634 } 00635 00636 private: 00637 /** 00638 * Initialize a lot of tensors. 00639 */ 00640 void setUp() 00641 { 00642 one( 0, 0 ) = one( 1, 1 ) = one( 2, 2 ) = 1.0; 00643 sone( 0, 0 ) = sone( 1, 1 ) = sone( 2, 2 ) = 1.0; 00644 00645 rdm1( 0, 0 ) = 2; 00646 rdm1( 0, 1 ) = 3; 00647 rdm1( 0, 2 ) = 1; 00648 rdm1( 1, 0 ) = 0; 00649 rdm1( 1, 1 ) = 4; 00650 rdm1( 1, 2 ) = 0; 00651 rdm1( 2, 0 ) = 5; 00652 rdm1( 2, 1 ) = -3; 00653 rdm1( 2, 2 ) = -7; 00654 00655 rdm2( 0, 0 ) = -1; 00656 rdm2( 0, 1 ) = -1; 00657 rdm2( 0, 2 ) = -1; 00658 rdm2( 1, 0 ) = 3; 00659 rdm2( 1, 1 ) = 0; 00660 rdm2( 1, 2 ) = -2; 00661 rdm2( 2, 0 ) = -1; 00662 rdm2( 2, 1 ) = -2; 00663 rdm2( 2, 2 ) = -3; 00664 00665 res1( 0, 0 ) = 6; 00666 res1( 0, 1 ) = -4; 00667 res1( 0, 2 ) = -11; 00668 res1( 1, 0 ) = 12; 00669 res1( 1, 1 ) = 0; 00670 res1( 1, 2 ) = -8; 00671 res1( 2, 0 ) = -7; 00672 res1( 2, 1 ) = 9; 00673 res1( 2, 2 ) = 22; 00674 00675 srdm1( 0, 0 ) = 2; 00676 srdm1( 0, 1 ) = 3; 00677 srdm1( 0, 2 ) = 1; 00678 srdm1( 1, 1 ) = -1; 00679 srdm1( 1, 2 ) = 4; 00680 srdm1( 2, 2 ) = -3; 00681 00682 srdm2( 0, 0 ) = -2; 00683 srdm2( 0, 1 ) = 5; 00684 srdm2( 0, 2 ) = 2; 00685 srdm2( 1, 1 ) = 0; 00686 srdm2( 1, 2 ) = -3; 00687 srdm2( 2, 2 ) = 1; 00688 00689 // copy symmetric tensor srdm1 to an asymmetric tensor 00690 res2 = srdm1; 00691 00692 res3( 0, 0 ) = 13; 00693 res3( 0, 1 ) = 7; 00694 res3( 0, 2 ) = -4; 00695 res3( 1, 0 ) = -3; 00696 res3( 1, 1 ) = 3; 00697 res3( 1, 2 ) = 13; 00698 res3( 2, 0 ) = 12; 00699 res3( 2, 1 ) = 14; 00700 res3( 2, 2 ) = -13; 00701 00702 res4( 0, 0 ) = 9; 00703 res4( 0, 1 ) = 15; 00704 res4( 0, 2 ) = -5; 00705 res4( 1, 0 ) = 26; 00706 res4( 1, 1 ) = -7; 00707 res4( 1, 2 ) = -25; 00708 res4( 2, 0 ) = -13; 00709 res4( 2, 1 ) = 28; 00710 res4( 2, 2 ) = 22; 00711 } 00712 00713 /** 00714 * A helper function that implements the simple approach to tensor evaluation. 00715 * 00716 * \param t The tensor. 00717 * \param v The gradient. 00718 * 00719 * \return value 00720 */ 00721 double calcTens( WTensorSym< 4, 3, double > const& t, WVector3d const& v ) 00722 { 00723 double res = 0.0; 00724 for( int a = 0; a < 3; ++a ) 00725 { 00726 for( int b = 0; b < 3; ++b ) 00727 { 00728 for( int c = 0; c < 3; ++c ) 00729 { 00730 for( int d = 0; d < 3; ++d ) 00731 { 00732 res += v[ a ] * v[ b ] * v[ c ] * v[ d ] * t( a, b, c, d ); 00733 } 00734 } 00735 } 00736 } 00737 return res; 00738 } 00739 00740 //! a test tensor 00741 WTensor< 2, 3, int > one; 00742 //! a test tensor 00743 WTensor< 2, 3, int > zero; 00744 //! a test tensor 00745 WTensor< 2, 3, int > rdm1; 00746 //! a test tensor 00747 WTensor< 2, 3, int > rdm2; 00748 //! a test tensor 00749 WTensor< 2, 3, int > res1; 00750 //! a test tensor 00751 WTensor< 2, 3, int > res2; 00752 //! a test tensor 00753 WTensor< 2, 3, int > res3; 00754 //! a test tensor 00755 WTensor< 2, 3, int > res4; 00756 //! a test tensor 00757 WTensorSym< 2, 3, int > sone; 00758 //! a test tensor 00759 WTensorSym< 2, 3, int > szero; 00760 //! a test tensor 00761 WTensorSym< 2, 3, int > srdm1; 00762 //! a test tensor 00763 WTensorSym< 2, 3, int > srdm2; 00764 }; 00765 00766 #endif // WTENSORFUNCTIONS_TEST_H