00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054 #ifdef UNITTEST
00055 #include "unit_test.h"
00056 using namespace matrix;
00057
00058 const D EPS=1e-9;
00059 bool comparetoidentity(const Matrix& m) {
00060 int worstdiagonal = 0;
00061 D maxunitydeviation = 0.0;
00062 D currentunitydeviation;
00063 for ( int i = 0; i < m.getM(); i++ ) {
00064 currentunitydeviation = m.val(i,i) - 1.;
00065 if ( currentunitydeviation < 0.0) currentunitydeviation *= -1.;
00066 if ( currentunitydeviation > maxunitydeviation ) {
00067 maxunitydeviation = currentunitydeviation;
00068 worstdiagonal = i;
00069 }
00070 }
00071 int worstoffdiagonalrow = 0;
00072 int worstoffdiagonalcolumn = 0;
00073 D maxzerodeviation = 0.0;
00074 D currentzerodeviation ;
00075 for ( int i = 0; i < m.getM(); i++ ) {
00076 for ( int j = 0; j < m.getN(); j++ ) {
00077 if ( i == j ) continue;
00078 currentzerodeviation = m.val(i,j);
00079 if ( currentzerodeviation < 0.0) currentzerodeviation *= -1.0;
00080 if ( currentzerodeviation > maxzerodeviation ) {
00081 maxzerodeviation = currentzerodeviation;
00082 worstoffdiagonalrow = i;
00083 worstoffdiagonalcolumn = j;
00084 }
00085
00086 }
00087 }
00088
00089
00090
00091
00092
00093 return (maxunitydeviation < EPS && maxzerodeviation < EPS);
00094 }
00095
00096
00097 UNIT_TEST_DEFINES
00098
00099 DEFINE_TEST( check_creation ) {
00100 cout << "\n -[ Creation and Access ]-\n";
00101 Matrix M1(3,3);
00102 D testdata[9]={1,0,0, 0,1,0, 0,0,1};
00103 Matrix M2(3,3);
00104 M2.set(testdata);
00105 M1.toId();
00106 unit_assert( "id=identity_set", M1 == M2 );
00107 D testdata2[12]={1,2,3, 4,5,6, 7,8,9, 10,11,12};
00108 Matrix M3(4,3, testdata2);
00109 unit_assert( "dimension of matrix", M3.getM() == 4 && M3.getN() == 3 );
00110 unit_assert( "field check",
00111 M3.val(0,2) == 3 && M3.val(2,1) == 8 && M3.val(3,0) == 10 );
00112 Matrix M4(M3);
00113 unit_assert( "copy constructor", M3 == M4 );
00114 Matrix M5(1,3,testdata2+3);
00115 unit_assert( "row", M3.row(1) == M5 );
00116
00117 unit_pass();
00118 }
00119
00120 DEFINE_TEST( check_vector_operation ) {
00121 cout << "\n -[ Vector Operations ]-\n";
00122 D testdata[3]={-1,3,2};
00123 const Matrix V1(1,3, testdata);
00124 const Matrix V2(3,1, testdata);
00125 Matrix V3(V1);
00126
00127 V3.toTranspose();
00128 unit_assert( "transpose", V3 == V2 );
00129 V3.toTranspose();
00130 unit_assert( "double transpose", V3 == V1 );
00131 D testdata2[3]={-2,6,4};
00132 Matrix V4(1,3,testdata2);
00133 V3.add(V1,V1);
00134 unit_assert( "add", V3 == V4 );
00135 D testdata3[3]={1,-3,-2};
00136 Matrix V5(1,3,testdata3);
00137 V4.sub(V1,V3);
00138 unit_assert( "sub", V4 == V5 );
00139
00140 D testdata4[3]={3,-9,-6};
00141 V5.set(testdata4);
00142 V4.copy(V1);
00143 V4.toMult(-3.0);
00144 unit_assert( "mult with scalar I", V4 == V5 );
00145 V4.mult(V1,-3.0);
00146 unit_assert( "mult with scalar II", V4 == V5 );
00147
00148 double f = 14;
00149 Matrix V6(1,1,&f);
00150 V3.copy(V1);
00151 V3.toTranspose();
00152 V4.mult(V1,V3);
00153 unit_assert( "scalarproduct", V4 == V6 );
00154 V4.copy(V1);
00155 unit_assert( "scalarproduct with exp(2)", V4.multMT() == V6 );
00156 V4.mult(V3,V1);
00157 unit_assert( "vector^T*vector=matrix", V4.getM() == 3);
00158 unit_assert( "vector^T*vector=exp(2)", V3.multMT().getM() == 3);
00159
00160 unit_pass();
00161 }
00162
00163 DEFINE_TEST( check_matrix_operation ) {
00164 cout << "\n -[ Matrix Operations ]-\n";
00165 D testdata[6]={1,2,3, 4,5,6 };
00166 const Matrix M1(2,3, testdata);
00167 D testdata2[6]={1,4, 2,5, 3,6 };
00168 const Matrix M2(3,2, testdata2);
00169 Matrix M3(M1);
00170
00171 M3.toTranspose();
00172 unit_assert( "transpose", M3 == M2 );
00173 D testdata3[6]={2,4,6, 8,10,12};
00174 Matrix M4(2,3,testdata3);
00175 M3.add(M1,M1);
00176 unit_assert( "add", M3 == M4 );
00177
00178 D testdata4[6]={-0.5, -1, -1.5, -2, -2.5, -3};
00179 M3.set(testdata4);
00180 M4.copy(M1);
00181 M4.toMult(-0.5);
00182 unit_assert( "mult with scalar", M4 == M3 );
00183
00184 D testdata5[6] = {2, -1, 0, 3, 2, -2};
00185 D testdata6[4] = {8, -1, 20, -1};
00186 Matrix M5 (3,2, testdata5);
00187 Matrix M6(2,2,testdata6);
00188 M4.mult(M1,M5);
00189 unit_assert( "mult(matrix, matrix)", M4 == M6 );
00190
00191 M3.copy(M1);
00192 M4.copy(M1);
00193 M4.toExp(1);
00194 unit_assert( "exp(1)", M3 == M4 );
00195 M3.toTranspose();
00196 M4.toExp(T);
00197 unit_assert( "exp(T)=transpose", M3 == M4 );
00198 M3.toId();
00199 M4.toExp(0);
00200 unit_assert( "exp(0)=id", M3 == M4 );
00201
00202 D testdata7[16] = {1,2,3,4, -4,2,1,3, 0.3,-0.9, 4, -3, 1,0.5,0.3,5.0};
00203 Matrix M7(4,4,testdata7);
00204 Matrix M8(M7);
00205 M7.toExp(-1);
00206 M4.mult(M7,M8);
00207 unit_assert( "exp(-1)*exp(1)=id", comparetoidentity(M4) );
00208
00209 D testdata9[6] = {sin(1.0),sin(2.0),sin(3.0), sin(4.0),sin(5.0),sin(6.0) };
00210 Matrix M9(2,3,testdata9);
00211 M4.copy(M1);
00212 M4.toMap(sin);
00213 unit_assert( "map(sin)", M4 == M9 );
00214
00215 D testdata10[6] = {2,4,6, -0.4,-0.5,-0.6 };
00216 D testdata11[2] = {2,-0.1};
00217 Matrix M10(2,3,testdata10);
00218 Matrix M11(2,1,testdata11);
00219 M4.copy(M1);
00220 M4.toMultrowwise(M11);
00221 unit_assert( "multrowwise()", M4 == M10 );
00222 D testdata12[6] = {2,1,0, 8, 2.5, 0 };
00223 D testdata13[3] = {2, 0.5, 0};
00224 M10.set(2,3,testdata12);
00225 Matrix M12(3,1,testdata13);
00226 M4.copy(M1);
00227 M4.toMultcolwise(M12);
00228 unit_assert( "multrowwise()", M4 == M10 );
00229
00230 M3.copy(M1);
00231 M4.copy(M1);
00232 M4.toTranspose();
00233 M5 = M3.multMT();
00234 M6.mult(M1,M4);
00235 unit_assert( "multMT() ", M5 == M6 );
00236 M5 = M3.multTM();
00237 M6.mult(M4,M1);
00238 unit_assert( "multTM() ", M5 == M6 );
00239
00240 D testdata20[12]={1,2,3, 4,5,6, 1,2,3, 4,5,6 };
00241 const Matrix M20(4,3, testdata20);
00242 const Matrix M21 = M1.above(M1);
00243 unit_assert( "above() ", M20 == M21 );
00244
00245
00246
00247 unit_pass();
00248 }
00249
00250 DEFINE_TEST( check_matrix_operators ) {
00251 cout << "\n -[ Matrix Operators (+ - * ^)]-\n";
00252 D testdata[6]={1,2,3, 4,5,6 };
00253 const Matrix M1(2,3, testdata);
00254 D testdata2[6]={1,4, 2,5, 3,6 };
00255 const Matrix M2(3,2, testdata2);
00256 unit_assert( "^T ", (M1^T) == M2 );
00257 D testdata3[6]={2,4,6, 8,10,12};
00258 Matrix M4(2,3,testdata3);
00259 unit_assert( "+ ", M1+M1 == M4 );
00260 unit_assert( "- ", M1+M1-M1 == M1 );
00261
00262 D testdata4[6]={-0.5, -1, -1.5, -2, -2.5, -3};
00263 Matrix M3(2,3, testdata4);
00264 unit_assert( "* scalar", M1*(-0.5) == M3 );
00265
00266 D testdata5[6] = {2, -1, 0, 3, 2, -2};
00267 D testdata6[4] = {8, -1, 20, -1};
00268 Matrix M5 (3,2, testdata5);
00269 Matrix M6(2,2,testdata6);
00270 unit_assert( "* ", M1*M5 == M6 );
00271
00272 unit_assert( "^1 ", (M1^1) == M1 );
00273 M3.toId();
00274 unit_assert( "^0=id ", (M1^0) == M3 );
00275
00276 D testdata7[16] = {1,2,3,4, -4,2,1,3, 0.3,-0.9, 4, -3, 1,0.5,0.3,5.0};
00277 Matrix M7(4,4,testdata7);
00278 unit_assert( "^1 * ^-1=id ", comparetoidentity(M7*(M7^-1)) );
00279 unit_pass();
00280 }
00281
00282 DEFINE_TEST( speed ) {
00283 cout << "\n -[ Speed: Inverion]-\n";
00284 #ifndef NDEBUG
00285 cout << " DEBUG MODE! use -DNDEBUG -O3 (not -g) to get full performance\n";
00286 #endif
00287 Matrix M1;
00288 srand(time(0));
00289 D testdata0[9] = {1,2, -4,2};
00290 Matrix M2(2,2,testdata0);
00291 UNIT_MEASURE_START("2x2 Matrix inversion", 100000)
00292 M1 = (M2^-1);
00293 UNIT_MEASURE_STOP("");
00294 unit_assert( "validation", comparetoidentity(M1*M2));
00295
00296
00297 D testdata1[9] = {1,2,3, -4,2,1, 0.3,-0.9};
00298 Matrix M3(3,3,testdata1);
00299 UNIT_MEASURE_START("3x3 Matrix inversion", 100000)
00300 M1 = (M3^-1);
00301 UNIT_MEASURE_STOP("");
00302 unit_assert( "validation", comparetoidentity(M1*M3));
00303
00304 D testdata2[16] = {1,2,3,4, -4,2,1,3, 0.3,-0.9, 4, -3, 1,0.5,0.3,5.0};
00305 Matrix M4(4,4,testdata2);
00306 UNIT_MEASURE_START("4x4 Matrix inversion", 100000)
00307 M1 = (M4^-1);
00308 UNIT_MEASURE_STOP("");
00309 unit_assert( "validation", comparetoidentity(M1*M4));
00310
00311 Matrix M20(20,20);
00312 for (int i=0; i < M20.getM(); i++)
00313 for (int j=0; j < M20.getN(); j++) {
00314 M20.val(i,j) = -22+(100. * rand())/RAND_MAX;
00315 }
00316 UNIT_MEASURE_START("20x20 Matrix inversion",1000)
00317 M1 = (M20^-1);
00318 UNIT_MEASURE_STOP("");
00319 unit_assert( "validation", comparetoidentity(M1*M20));
00320
00321 Matrix M200(200,200);
00322 rand();
00323 for (int i=0; i < M200.getM(); i++)
00324 for (int j=0; j < M200.getN(); j++) {
00325 M200.val(i,j) = -22+(100. * rand())/RAND_MAX;
00326 }
00327 UNIT_MEASURE_START("200x200 Matrix inversion",2)
00328 M1 = (M200^-1);
00329 UNIT_MEASURE_STOP("");
00330 unit_assert( "validation", comparetoidentity(M1*M200));
00331
00332 cout << "\n -[ Speed: Other Operations]-\n";
00333 UNIT_MEASURE_START("20x20 Matrix multiplication with assignment",5000)
00334 M1 = M20*M20;
00335 UNIT_MEASURE_STOP("");
00336 UNIT_MEASURE_START("20x20 Matrix addition with assignment",100000)
00337 M1= (M1 + M20);
00338 UNIT_MEASURE_STOP("");
00339 UNIT_MEASURE_START("20x20 Matrix inplace addition",100000)
00340 M1 += M20;
00341 UNIT_MEASURE_STOP("");
00342 UNIT_MEASURE_START("20x20 Matrix transposition",100000)
00343 M1 += M20;
00344 UNIT_MEASURE_STOP("");
00345
00346 unit_pass();
00347 }
00348
00349 UNIT_TEST_RUN( "Matrix Tests" )
00350 ADD_TEST( check_creation )
00351 ADD_TEST( check_vector_operation )
00352 ADD_TEST( check_matrix_operation )
00353 ADD_TEST( check_matrix_operators )
00354 ADD_TEST( speed )
00355
00356 UNIT_TEST_END
00357
00358 #endif
