00001 /*
00002 File: VLTest.cc
00003
00004 Function: Test program for the VL class.
00005
00006 Author(s): Andrew Willmott
00007
00008 Copyright: (c) 1995-2000, Andrew Willmott
00009
00010 Notes:
00011 */
00012
00013 #include <stdio.h>
00014 #include <stdlib.h>
00015 #include <iostream.h>
00016 #include <iomanip.h>
00017
00018 #include "vl/VLTest.h"
00019
00020 // --- Prototypes -------------------------------------------------------------
00021
00022 Void TestBasics();
00023
00024 Void Test2DStuff();
00025 Void Test3DStuff();
00026 Void Test4DStuff();
00027
00028 Void TestHStuff2D();
00029 Void TestHStuff3D();
00030
00031 Void TestND();
00032 Void TestNDSub();
00033 Void TestNDNumerical();
00034
00035 Void TestSparse();
00036 Void TestSparseSub();
00037 Void TestSparseNumerical();
00038
00039
00040 Void TestInit();
00041 Void ComplexTest();
00042
00043 #define TYPE_NAME(X) type_name((X*)0)
00044 inline char *type_name(Float*) { return("Float"); };
00045 inline char *type_name(Double*) { return("Double"); };
00046 inline char *type_name(Int*) { return("Int"); };
00047
00048
00049 int main(int, char **)
00050 {
00051 cout << "Testing VL library, version " << VL_VERSION << endl;
00052 cout << "Real type: " << TYPE_NAME(Real) << endl;
00053
00054 cout << "----------------------------------" << endl;
00055
00056 Test2DStuff();
00057 Test3DStuff();
00058 Test4DStuff();
00059
00060 TestHStuff2D();
00061 TestHStuff3D();
00062
00063 TestND();
00064 TestNDSub();
00065 TestNDNumerical();
00066
00067 TestSparse();
00068 TestSparseSub();
00069 TestSparseNumerical();
00070
00071 #ifdef VL_COMPLEX
00072 ComplexTest();
00073 #endif
00074 TestInit();
00075
00076 cout << "\n\n--- Finished! ---" << endl;
00077
00078 return(0);
00079 }
00080
00081
00082 // --- Test routines ----------------------------------------------------------
00083
00084 Void TestHStuff2D()
00085 {
00086 Vec2r x;
00087 Vec3r y;
00088 Mat2s M;
00089
00090 cout << "\n+ TestHStuff2D\n\n";
00091
00092 x = Vec2r(1,2);
00093 cout << "x is: " << x << endl;
00094 M = Rot2s(vl_halfPi);
00095 cout << "rot(pi/2) is: " << Rot2s(vl_halfPi) << endl;
00096 x = Rot2s(vl_halfPi) * x;
00097 cout << "x after rot(pi/2) is: " << x << endl;
00098 x = Scale2s(Vec2s(0.3, 0.2)) * x;
00099 cout << "x after scale(0.3, 0.2) is: " << x << endl;
00100
00101 y = Vec3r(x, 0.5);
00102 cout << "y is: " << y << endl;
00103 x = proj(y);
00104 cout << "proj(y) is: " << x << endl;
00105
00106 x = proj(HRot3s(1.3) * HTrans3s(Vec2s(1,1)) *
00107 HScale3s(Vec2s(1,2)) * Vec3r(x, 1));
00108 cout << "HRot3s(1.3) * HTrans3s(Vec2s(1,1)) * HScale3s(Vec2s(1,2)) * y = "
00109 << x << endl;
00110 }
00111
00112 Void TestHStuff3D()
00113 {
00114 Vec3r x;
00115 Vec4r y;
00116 Mat4s z;
00117
00118 cout << "\n+ TestHStuff3D\n\n";
00119
00120 x = Vec3r(1,2,3);
00121
00122 cout << "rot(pi/2, vl_x) is: " << Rot3s(vl_x, vl_halfPi) << endl;
00123 x = x * Rot3s(vl_x, vl_halfPi);
00124 cout << "x after rot(pi/2, vl_x) is: " << x << endl;
00125 x = x * Scale3s(Vec3s(0.3, 0.2, 0.3));
00126 cout << "x after scale(0.3, 0.2, 0.3) is: " << x << endl;
00127
00128 y = Vec4r(x, 0.5);
00129 cout << "y is: " << y << endl;
00130 x = proj(y);
00131 cout << "proj(y) is: " << x << endl;
00132
00133 x = proj(HRot4s(vl_x, 1.3) * HTrans4s(vl_1) * HScale4s(Vec3s(1,2,1)) * y);
00134 cout << "HRot4s(vl_x, 1.3) * HTrans4s(vl_1) "
00135 "* HScale4s(Vec3s(1,2,1)) * y = " << x;
00136 }
00137
00138
00139 Void Test2DStuff()
00140 {
00141 Vec2r x(1,2);
00142 Vec2r y(5,6);
00143
00144 cout << "\n+ Test2DStuff\n\n";
00145
00146 cout << "x: " << x << ", y: " << y << "\n\n";
00147
00148 cout << "x + y * (y * x * 2) : " << x + y * (y * x * 2) << endl;
00149 cout << "x dot y : " << dot(x, y) << endl;
00150
00151 cout << "cross(x) : " << cross(x) << endl;
00152 cout << "len : " << len(x) << endl;
00153 cout << "sqrlen : " << sqrlen(x) << endl;
00154 cout << "norm : " << norm(x) << endl;
00155 cout << "len of norm : " << len(norm(x)) << "\n\n";
00156
00157 Mat2s M(1,2,3,4);
00158 Mat2s N; N.MakeDiag(2.0);
00159 Mat2s I = vl_I;
00160
00161 cout << "M : " << M << endl;
00162
00163 cout << "M * x : " << M * x << endl;
00164 cout << "x * M : " << x * M << endl;
00165
00166 cout << "adj : " << adj(M) << endl;
00167 cout << "det : " << det(M) << endl;
00168 cout << "trace : " << trace(M) << endl;
00169 cout << "inv : \n" << inv(M) << endl;
00170 cout << "M * inv : \n" << clamped(M * inv(M)) << "\n" << endl;
00171
00172 cout << "Vec2 consts: " << Vec2r(vl_0) << Vec2r(vl_x)
00173 << Vec2r(vl_y) << Vec2r(vl_1) << endl;
00174 cout << "Mat2 consts:\n" << Mat2s(vl_Z) << endl << Mat2s(vl_I)
00175 << endl << Mat2s(vl_B) << "\n\n";
00176
00177 M = Rot2s(1.3) * Scale2s(Vec2s(2,1));
00178
00179 cout << "M : \n" << M << endl;
00180
00181 cout << "M * x : " << M * x << endl;
00182 cout << "x * M : " << x * M << endl;
00183
00184 cout << "adj : " << adj(M) << endl;
00185 cout << "det : " << det(M) << endl;
00186 cout << "trace : " << trace(M) << endl;
00187 cout << "inv : \n" << inv(M) << endl;
00188 cout << "M * inv : \n" << clamped(M * inv(M)) << "\n" << endl;
00189 }
00190
00191 Void Test3DStuff()
00192 {
00193 Vec3r x(1,2,3);
00194 Vec3r y(5,6,7);
00195
00196 cout << "\n+ Test3DStuff\n\n";
00197
00198 cout << "x: " << x << ", y: " << y << "\n\n";
00199
00200 cout << "x + y * (y * x * 2) : " << x + y * (y * x * 2) << endl;
00201 cout << "x dot y : " << dot(x, y) << endl;
00202
00203 cout << "cross(x,y) : " << cross(x,y) << endl;
00204 cout << "cross(x, y) . x : " << dot(cross(x, y), x) << endl;
00205 cout << "cross(x, y) . y : " << dot(cross(x, y), y) << endl;
00206 cout << "len : " << len(x) << endl;
00207 cout << "sqrlen : " << sqrlen(x) << endl;
00208 cout << "norm : " << norm(x) << endl;
00209 cout << "len of norm : " << len(norm(x)) << "\n\n";
00210
00211 Mat3s M(1,2,3,3,2,1,2,1,3);
00212 Mat3s N; N.MakeDiag(2.0);
00213 Mat3s I = vl_I;
00214
00215 cout << "M : \n" << M << endl;
00216
00217 cout << "M * x : " << M * x << endl;
00218 cout << "x * M : " << x * M << endl;
00219
00220 cout << "adj : " << adj(M) << endl;
00221 cout << "det : " << det(M) << endl;
00222 cout << "trace : " << trace(M) << endl;
00223 cout << "inv : \n" << inv(M) << endl;
00224 cout << "M * inv : \n" << clamped(M * inv(M)) << endl;
00225
00226 cout << "Vec3 consts: " << Vec3r(vl_0) << Vec3r(vl_x)
00227 << Vec3r(vl_y) << Vec3r(vl_z) << Vec3r(vl_1) << endl;
00228 cout << "Mat3 consts:\n" << Mat3s(vl_Z) << endl << Mat3s(vl_I)
00229 << endl << Mat3s(vl_B) << "\n\n";
00230
00231 M = Rot3s(vl_y, 1.3) * Scale3s(Vec3s(2,4,2));
00232
00233 cout << "M :\n" << M << endl;
00234
00235 cout << "M * x : " << M * x << endl;
00236 cout << "x * M : " << x * M << endl;
00237
00238 cout << "adj : " << adj(M) << endl;
00239 cout << "det : " << det(M) << endl;
00240 cout << "trace : " << trace(M) << endl;
00241 cout << "inv : \n" << inv(M) << endl;
00242 cout << "M * inv : \n" << clamped(M * inv(M)) << endl;
00243 }
00244
00245 Void Test4DStuff()
00246 {
00247 Vec4r x(1,2,3,4);
00248 Vec4r y(5,6,7,8);
00249 Vec4r z(1,0,0,0);
00250
00251 cout << "\n+ Test4DStuff\n\n";
00252
00253 cout << "x: " << x << ", y: " << y << ", z: " << z << "\n\n";
00254
00255 cout << "x + y * (z * x * 2) : " << x + y * (z * x * 2) << endl;
00256 cout << "x dot y : " << dot(x, y) << "\n\n";
00257
00258 cout << "cross(x,y,z) : " << cross(x,y,z) << endl;
00259 cout << "cross(x,y,z).x : " << dot(cross(x,y,z), x) << endl;
00260 cout << "cross(x,y,z).y : " << dot(cross(x,y,z), y) << endl;
00261 cout << "cross(x,y,z).z : " << dot(cross(x,y,z), z) << endl;
00262 cout << "len : " << len(x) << endl;
00263 cout << "sqrlen : " << sqrlen(x) << endl;
00264 cout << "norm : " << norm(x) << endl;
00265 cout << "len of norm : " << len(norm(x)) << "\n\n";
00266
00267
00268 Mat4s M(1,2,3,0, 2,3,0,5, 3,0,5,6, 0,5,6,7);
00269 Mat4s N; N.MakeBlock(2.0);
00270 Mat4s I = vl_I;
00271
00272 cout << "M : \n" << M << endl;
00273
00274 cout << "M * x : " << M * x << endl;
00275 cout << "x * M : " << x * M << endl;
00276
00277 cout << "adj : " << adj(M) << endl;
00278 cout << "det : " << det(M) << endl;
00279 cout << "trace : " << trace(M) << endl;
00280 cout << "inv : \n" << inv(M) << endl;
00281 cout << "M * inv : \n" << clamped(M * inv(M)) << endl;
00282
00283 cout << "Vec4 consts: " << Vec4r(vl_0) << Vec4r(vl_x) << Vec4r(vl_y)
00284 << Vec4r(vl_z) << Vec4r(vl_w) << Vec4r(vl_1) << endl;
00285 cout << "Mat4 consts:\n" << Mat4s(vl_Z) << endl << Mat4s(vl_I) << endl
00286 << Mat4s(vl_B) << "\n\n";
00287
00288 M = HScale4s(Vec3s(2,3,4));
00289 M *= HRot4s(vl_y, 1.256);
00290
00291 cout << "M : " << M << endl;
00292
00293 cout << "M * x : " << M * x << endl;
00294 cout << "x * M : " << x * M << endl;
00295
00296 cout << "adj : " << adj(M) << endl;
00297 cout << "det : " << det(M) << endl;
00298 cout << "trace : " << trace(M) << endl;
00299 cout << "inv : \n" << inv(M) << endl;
00300 cout << "M * inv : \n" << clamped(M * inv(M)) << endl;
00301
00302 }
00303
00304 ostream &operator << (ostream &s, const TVec &v);
00305
00306 Void TestND()
00307 {
00308 cout << "\n+ TestND\n\n";
00309
00310 Vecr x(4, 1.0, 2.0, 3.0, 4.0);
00311 Vecr sx(3, 1.0, 2.0, 3.0);
00312 Vecr y(4, 5.0, 6.0, 7.0, 8.0);
00313 Vecr z(4, 4.0, 3.0, 2.0, 1.0);
00314
00315 Mats M(4, 3, 4.0, 3.0, 2.0, 1.0, 4.0, 3.0, 2.0, 1.0, 5.0,
00316 2.0, 1.0, 4.0, 3.0, 2.0, 1.0, 5.0);
00317 Mats N(4, 3, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,
00318 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0);
00319
00320 cout << "M:\n" << M;
00321 cout << "\nN:\n" << N;
00322 cout << "\ntrans(N):\n" << trans(N);
00323 cout << "\nM + N:\n" << M + N;
00324 cout << "\nM * trans(N):\n" << M * trans(N);
00325 cout << "x: " << x << ", sx: " << sx << endl;
00326
00327 z = x + y;
00328 cout << "z: " << z << endl;
00329
00330 cout << "M * sx : " << (M * sx) << endl;
00331 cout << "x * M : " << x * M << endl;
00332
00333 cout << "x: " << x << ", y: " << y << ", z: " << z << endl;
00334
00335 cout << "x + y: " << x + y << endl;
00336 cout << "x * y: " << x * y << endl;
00337 cout << "x dot y: " << dot(x, y) << endl;
00338 cout << "x + (y dot z) * x * 2 : " << x + (dot(y, z) * x * 2.0) << endl;
00339 cout << "x + (y dot z) * x * 2 : " << ((float)2.0 * dot(y, z)) * x << endl;
00340 cout << "x + (y dot z) * x * 2 : " << x + dot(y, z) * x << endl;
00341
00342 cout << "len : " << len(x) << endl;
00343 cout << "sqrlen : " << sqrlen(x) << endl;
00344 cout << "norm : " << norm(x) << endl;
00345 cout << "len of norm : " << len(norm(x)) << endl;
00346 }
00347
00348 #ifdef OLD
00349 Void MatTest()
00350 {
00351 Matf m(50, 50);
00352 Int i, j;
00353
00354 for (i = 0; i < m.Rows(); i++)
00355 for (j = 0; j < m.Cols(); j++)
00356 m[i][j] = i * 100 + j;
00357
00358 cout << row(sub(sub(m, 20, 30, 10, 10), 3, 3, 4, 4), 2) << endl;
00359 row(sub(sub(m, 20, 30, 10, 10), 3, 3, 4, 4), 2) = Vecf(4, vl_1);
00360 cout << row(sub(sub(m, 20, 30, 10, 10), 3, 3, 4, 4), 2) << endl;
00361
00362 cout << col(sub(sub(m, 20, 30, 10, 10), 3, 3, 4, 4), 2) << endl;
00363 col(sub(sub(m, 20, 30, 10, 10), 3, 3, 4, 4), 2) = 2 * Vecf(4, vl_1);
00364 cout << col(sub(sub(m, 20, 30, 10, 10), 3, 3, 4, 4), 2) << endl;
00365
00366 cout << diag(sub(sub(m, 20, 30, 10, 10), 3, 3, 4, 4), 2) << endl;
00367 diag(sub(sub(m, 20, 30, 10, 10), 3, 3, 4, 4), 2) = Vecf(2, vl_1);
00368 cout << diag(sub(sub(m, 20, 30, 10, 10), 3, 3, 4, 4), 2) << endl;
00369
00370 cout << last(first(sub(m[0], 20, 30), 6), 2) << endl;
00371 last(first(sub(m[0], 20, 30), 6), 2) = Vecf(2, vl_1);
00372 cout << last(first(sub(m[0], 20, 30), 6), 2) << endl;
00373
00374 cout << m[40] << endl;
00375 ToWaveletTransformRow(row(m, 40));
00376 cout << m[40] << endl;
00377 }
00378 #endif
00379
00380 Void TestNDSub()
00381 {
00382 cout << "\n+ TestNDSub\n\n";
00383
00384 Vecr x(8, 1.0, 2.0, 3.0, 4.0, 1.0, 2.0, 3.0, 4.0);
00385 Vecr y(16, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0,
00386 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0);
00387 Vecr z;
00388
00389 cout << "x: " << x << ", y: " << y << endl;
00390
00391 cout << "sub(y, 3, 3): " << sub(y, 3, 3) << endl;
00392 z = sub(x, 2, 6) + sub(y, 10, 6);
00393 cout << "sub(x, 2, 6) + sub(y, 10, 6): " << z << endl;
00394
00395 sub(y, 5, 6) = Vecr(6, 88.0, 77.0, 66.0, 55.0, 44.0, 33.0);
00396 sub(x, 0, 2) = Vecr(2, 22.0, 11.0);
00397
00398 cout << "x: " << x << ", y: " << y << endl;
00399
00400 z = z + sub(y, 5, 6);
00401 sub(y, 5, 6) = sub(y, 5, 6) + z;
00402
00403 cout << "z: " << z << ", y: " << y << endl;
00404
00405 cout << "\n\n";
00406
00407 Mats M(10, 10, vl_I);
00408 Mats N(4, 3, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0,
00409 9.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0);
00410
00411 cout << "sub(N, 1, 1, 2, 2): " << endl << sub(N, 1, 1, 2, 2);
00412
00413 sub(M, 5, 3, 4, 3) = N;
00414
00415 cout << "\nM:\n" << M;
00416
00417 cout << "\nDiagonals of M: \n\n";
00418
00419 Int i;
00420
00421 for (i = 1 - M.Rows(); i < M.Cols(); i++)
00422 cout << diag(M, i) << endl;
00423
00424 cout << "\nCol 4 of M: \n" << col(M, 4);
00425
00426 diag(M, 0) = Vecs(10, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0);
00427 diag(M, 1) = Vecs(9, vl_1) * 2.0;
00428 diag(M, -1) = diag(M, 1) * 3.0;
00429
00430 cout << "\nM:\n" << M << endl;
00431 }
00432
00433 Void TestNDNumerical()
00434 {
00435 Mats P(4, 4,
00436 1.0, 2.0, 3.0, 0.0,
00437 2.0, 3.0, 0.0, 5.0,
00438 3.0, 0.0, 5.0, 6.0,
00439 0.0, 5.0, 6.0, 7.0
00440 );
00441 Mats Q;
00442 Mat4s P1(
00443 1.0, 2.0, 3.0, 0.0,
00444 2.0, 3.0, 0.0, 5.0,
00445 3.0, 0.0, 5.0, 6.0,
00446 0.0, 5.0, 6.0, 7.0
00447 );
00448 Mat4s Q1;
00449
00450 cout << "\n+ TestNDNumerical\n" << endl;
00451
00452 cout << "P:\n";
00453 cout << P;
00454
00455 cout << "\ninv(P):" << endl;
00456 Q = inv(P);
00457 cout << Q;
00458
00459 cout << "\nP * inv(P):\n";
00460 cout << clamped(P * Q);
00461
00462 cout << "\n\nReal:\n";
00463
00464 cout << "P1: " << P1 << endl;
00465 cout << "\ninv(P1):\n";
00466 Q1 = inv(P1);
00467 cout << Q1 << endl;
00468
00469 cout << "\nP1 * inv(P1): " << endl << clamped(P1 * Q1) << endl;
00470 cout << "\ninv(P) - inv(P1): " << endl << clamped(inv(P) - inv(P1));
00471 cout << endl << endl;
00472
00473 // --- SolveOverRelax -----------------------------------------------------
00474
00475 Mats A(4, 4,
00476 29.0, 2.0, 3.0, 0.0,
00477 2.0, 29.0, 0.0, 5.0,
00478 3.0, 0.0, 29.0, 6.0,
00479 0.0, 5.0, 6.0, 29.0);
00480
00481 Vecr x;
00482 Vecr b;
00483 Real error = 1.0;
00484 Int i = 0;
00485
00486 b.SetSize(4);
00487 b = A * Vecr(Vec4r(1.0, 2.0, 3.0, 4.0));
00488
00489 x = b;
00490 cout << "Solving Ax = b with over-relaxation..." << endl;
00491 cout << "A: \n" << A << endl;
00492 cout << "b: " << b << endl;
00493 cout << "start x: " << x << endl;
00494
00495 error = SolveOverRelax(A, x, b, 1e-8, 1.1);
00496
00497 // cout << "iterations: " << i << endl;
00498 cout << "x: " << x << endl;
00499 cout << "Ax - b: " << A * x - b << endl;
00500 cout << "Returned error: " << error << ", real error: " << sqrlen(A * x - b) << endl;
00501 cout << endl;
00502
00503 // --- SolveConjGrad ------------------------------------------------------
00504
00505 x = b;
00506
00507 cout << "Solving Ax = b with conjugate-gradient..." << endl;
00508 cout << "A:\n" << A;
00509 cout << "b: " << b << endl;
00510 cout << "start x: " << x << endl;
00511
00512 error = SolveConjGrad(A, x, b, 1e-8, &i);
00513
00514 cout << "iterations: " << i << endl;
00515 cout << "x: " << x << endl;
00516 cout << "Ax - b: " << A * x - b<< endl;
00517 cout << "Returned error: " << error << ", real error: " <<
00518 sqrlen(A * x - b) << endl;
00519 }
00520
00521
00522 Void TestSparse()
00523 {
00524 cout << "\n+ TestSparse\n" << endl;
00525
00526 SparseVecr v(30);
00527 SparseVecr z(30);
00528 SparseVecs w(Vecs(6, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0));
00529 Vecr x(30);
00530 Vecr y(30);
00531 Int i;
00532 Int i1[] = { 3, 4, 7, 19, VL_SV_END };
00533 Real e1[] = { 4.0, 5.0, 6.0, 7.0 };
00534
00535 for (i = 0; i < 30; i++)
00536 {
00537 x[i] = 1.0;
00538 y[i] = 2.0;
00539 }
00540
00541 v.Begin();
00542 v.AddElt(1, 3);
00543 v.AddElt(4, 2);
00544 v.AddElt(11, 3);
00545 v.AddElt(20, 2);
00546 v.End();
00547
00548 cout << v << endl;
00549 z = v;
00550
00551 v.SetElts(1, 4.0, 4, 0.0, 20, 3.0, 14, 6.0, VL_SV_END);
00552
00553 cout << v << endl;
00554
00555 v.Begin();
00556 for (i = 0; i < 100; i++)
00557 v.AddElt(i, i);
00558 v.End();
00559
00560 v.SetCompactPrint(true);
00561 cout << "compact: " << v << endl;
00562 v.SetCompactPrint(false);
00563
00564 cout << "v: " << v << endl;
00565 cout << "w: " << w << endl;
00566 cout << "x: " << x << endl;
00567 cout << "y: " << y << endl;
00568 cout << "z: " << z << endl;
00569
00570 v = 2 * v + x - z;
00571
00572 cout << "2 * v + x - z: " << v << endl;
00573 cout << "v * x: " << v * x << endl;
00574
00575 SparseMats M(20, 20, vl_I);
00576 SparseMats N(20, 20, vl_I);
00577 SVIters j;
00578
00579 w.SetSize(20);
00580 w.SetElts(5, 3.96, 10, 2.0, 15, 8.0, 19, 2.0, VL_SV_END);
00581
00582 M[10] = M[10].Overlay(w);
00583 M[15].Set(5, 6.0);
00584
00585 cout << M[15].Get(5) << endl;
00586
00587 cout << SparseVecr(20, 5, 15.0, 9, 21.0, 14, 20.0, VL_SV_END) << endl;
00588 cout << SparseVecs(20, i1, e1) << endl;
00589
00590 M[8].Overlay(SparseVecs(20, 5, 15.0, 9, 21.0, 14, 20.0, VL_SV_END));
00591 cout << M[8] << endl;
00592
00593 cout << "M:\n" << M;
00594
00595 N = M * M;
00596
00597 cout << "N = M^2:\n" << N;
00598 cout << "N - M:\n" << N - M;
00599
00600 }
00601
00602 Void TestSparseSub()
00603 {
00604 cout << "\n+ TestSparseSub\n\n";
00605
00606 SparseVecr x(Vecr(8, 1.0, 2.0, 3.0, 4.0, 1.0, 2.0, 3.0, 4.0));
00607 SparseVecr y(Vecr(16, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0,
00608 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0));
00609 SparseVecr z;
00610
00611 cout << "x: " << x << ", y: " << y << endl;
00612
00613 cout << "sub(y, 3, 3): " << sub(y, 3, 3) << endl;
00614 z = sub(x, 2, 6) + sub(y, 10, 6);
00615 cout << "z = sub(x, 2, 6) + sub(y, 10, 6): " << z << endl;
00616
00617 sub(y, 5, 6) = Vecr(6, 88.0, 77.0, 66.0, 55.0, 44.0, 33.0);
00618 sub(x, 0, 2) = Vecr(2, 22.0, 11.0);
00619
00620 cout << "x: " << x << ", y: " << y << endl;
00621
00622 z = z + sub(y, 5, 6);
00623 sub(y, 5, 6) = sub(y, 5, 6) + z;
00624
00625 cout << "z: " << z << ", y: " << y << endl;
00626
00627 cout << "\n\n";
00628
00629 SparseMats M(10, 10, vl_I);
00630 SparseMats N(Mats(4, 3, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0,
00631 9.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0));
00632
00633 cout << "sub(N, 1, 1, 2, 2): " << endl << sub(N, 1, 1, 2, 2);
00634
00635 sub(M, 5, 3, 4, 3) = N;
00636
00637 cout << "\nM:\n" << M;
00638
00639 cout << "\nDiagonals of M: \n\n";
00640
00641 Int i;
00642
00643 for (i = 1 - M.Rows(); i < M.Cols(); i++)
00644 cout << diag(M, i) << endl;
00645
00646 cout << "\nCol 4 of M: \n" << col(M, 4);
00647
00648 diag(M, 0) = Vecs(10, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0);
00649 diag(M, 1) = Vecs(9, vl_1) * 2.0;
00650 diag(M, -1) = diag(M, 1) * 3.0;
00651
00652 cout << "\nM:\n" << M << endl;
00653 }
00654
00655 Void TestSparseNumerical()
00656 {
00657 #ifdef VL_BROKEN
00658 XXX
00659 static Real sparseData[] =
00660 {
00661 0, 1, 1.0, 2, 2.0, 3, 4.0, VL_SV_END,
00662 1, 1, 2.0, 2, 3.0, 4, 4.0, VL_SV_END,
00663 2, 1, 3.0, 2, 4.0, 3, 6.0, VL_SV_END,
00664 3, 1, 5.0, 2, 6.0, 3, 7.0, VL_SV_END,
00665 VL_SV_END
00666 }
00667 static Real sparseData[] =
00668 {
00669 0, 1, 1.0, 2, 2.0, 3, 4.0, VL_SV_END,
00670 1, 1, 2.0, 2, 3.0, 4, 4.0, VL_SV_END,
00671 2, 1, 3.0, 2, 4.0, 3, 6.0, VL_SV_END,
00672 3, 1, 5.0, 2, 6.0, 3, 7.0, VL_SV_END,
00673 VL_SV_END
00674 }
00675 #endif
00676
00677 SparseMats P = Mats(4, 4, 1.0, 2.0, 3.0, 0.0,
00678 2.0, 3.0, 0.0, 5.0,
00679 3.0, 0.0, 5.0, 6.0,
00680 0.0, 5.0, 6.0, 7.0);
00681 SparseMats Q;
00682
00683
00684 cout << "\n+ TestSparseNumerical\n" << endl;
00685
00686 cout << "P:\n";
00687 cout << P;
00688
00689 cout << "\ninv(P):" << endl;
00690 Q = inv(P);
00691 cout << Q;
00692
00693 cout << "\nP * inv(P):\n";
00694 cout << P * Q;
00695
00696 cout << "\n\nReal:\n";
00697
00698 Mat4s P1(1,2,3,0, 2,3,0,5, 3,0,5,6, 0,5,6,7);
00699 Mats Q1;
00700
00701 cout << "P1: " << P1 << endl;
00702
00703 cout << "\ninv(P1):\n";
00704 Q1 = inv(Mats(P1));
00705 cout << Q1 << endl;
00706
00707 cout << "\nP1 * inv(P1): " << endl << clamped(P1 * Q1) << endl;
00708
00709
00710 cout << "\ninv(P) - inv(P1): " << endl << setprecision(4)
00711 << Q - SparseMats(Q1);
00712
00713 cout << endl <<
00714 "---------------------------------------------------------" << endl;
00715
00716 // --- SolveOverRelax -----------------------------------------------------
00717
00718 SparseMats A = Mats(4, 4,
00719 29.0, 2.0, 3.0, 0.0,
00720 2.0, 29.0, 0.0, 5.0,
00721 3.0, 0.0, 29.0, 6.0,
00722 0.0, 5.0, 6.0, 29.0);
00723
00724 Vecr x, b = A * Vecr(4, 1.0, 2.0, 3.0, 4.0);
00725 Real error = 1;
00726 Int i = 0;
00727
00728 x = b;
00729
00730 cout << "Solving Ax = b with over-relaxation..." << endl;
00731 cout << "A: \n" << A << endl;
00732 cout << "b: " << b << endl;
00733
00734 error = SolveOverRelax(A, x, b, 1e-8, 1.1);
00735
00736 cout << "x: " << x << endl;
00737 cout << "Ax: " << A * x << endl;
00738 cout << "Ax - b: " << A * x - b << endl;
00739 cout << "Returned error: " << error << ", real error: "
00740 << sqrlen(A * x - b) << endl;
00741 cout << endl;
00742
00743 // --- SolveConjGrad ------------------------------------------------------
00744
00745 x = b;
00746
00747 cout << "A:\n" << A;
00748 cout << "b: " << b << endl;
00749 cout << "x: " << x << endl;
00750
00751 cout << "Solving Ax = b with conjugate-gradient..." << endl;
00752
00753 error = SolveConjGrad(A, x, b, 1e-8, &i);
00754
00755 cout << "iterations: " << i << endl;
00756 cout << "x: " << x << endl;
00757 cout << "Ax - b: " << A * x - b<< endl;
00758 cout << "Returned error: " << error << ", real error: "
00759 << sqrlen(A * x - b) << endl;
00760 }
00761
00762
00763 Void TestBasics()
00764 {
00765 cout << "\n+ TestBasics\n\n";
00766
00767 Assert(false, "Bogus assertion");
00768 Expect(false, "This is true");
00769
00770 cout << "max is: " << Max(2.0, 5.0) << endl;
00771 cout << "min is: " << Min((Int) 2, (Int) 5) << endl;
00772 }
00773
00774 Void TestInit()
00775 {
00776 Vecr v00(10, vl_zero);
00777 Vecr v01(10, vl_one);
00778 Vecr v02(10, vl_x);
00779 Vecr v03(10, vl_axis(5));
00780 Vecr v04(10); v04.MakeBlock(5.0);
00781
00782 Mats m00(5, 5, vl_zero);
00783 Mats m01(5, 5, vl_one);
00784 Mats m02(5, 5); m02.MakeDiag(4.0);
00785 Mats m03(5, 5); m03.MakeBlock(2.2);
00786 Mats m04(5, 5); m04.MakeBlock(8.8);
00787
00788 SparseVecr sv00(10, vl_zero);
00789 // SparseVecr sv01(10, vl_one); illegal (would be stupid)
00790 SparseVecr sv02(10, vl_x);
00791 SparseVecr sv03(10, vl_axis(5));
00792
00793 SparseMats sm00(5, 5, vl_zero);
00794 SparseMats sm01(5, 5, vl_one);
00795 SparseMats sm02(5, 5); sm02.MakeDiag(4.0);
00796 // SparseMats sm03(5, 5); sm03.MakeBlock(2.2); illegal (ditto)
00797
00798 sub(m04, 2, 2, 2, 2) = Mat2s(vl_B) * 111.0;
00799
00800 cout << "--- Test init code -------------------------------------------"
00801 << endl << endl;
00802
00803 cout << v00 << endl;
00804 cout << v01 << endl;
00805 cout << v02 << endl;
00806 cout << v03 << endl;
00807 cout << v04 << endl;
00808
00809 cout << endl;
00810
00811 cout << m00 << endl;
00812 cout << m01 << endl;
00813 cout << m02 << endl;
00814 cout << m03 << endl;
00815 cout << m04 << endl;
00816 cout << endl;
00817
00818 cout << sv00 << endl;
00819 cout << sv02 << endl;
00820 cout << sv03 << endl;
00821
00822 cout << endl;
00823
00824 cout << sm00 << endl;
00825 cout << sm01 << endl;
00826 cout << sm02 << endl;
00827 cout << endl;
00828
00829 cout << "testing comparisons..." << endl;
00830 cout << endl;
00831 cout << (Mat2s(vl_0) == vl_0) << endl;
00832 cout << (Mat3s(vl_0) == vl_0) << endl;
00833 cout << (Mat4s(vl_0) == vl_0) << endl;
00834 cout << (Mats(8, 8, vl_0) == Mats(8, 8, vl_0)) << endl;
00835
00836 cout << (Mat2s(vl_1) == vl_0) << endl;
00837 cout << (Mat3s(vl_1) == vl_0) << endl;
00838 cout << (Mat4s(vl_1) == vl_0) << endl;
00839 cout << (Mats(8, 8, vl_1) == Mats(8, 8, vl_0)) << endl;
00840 }
00841
00842
00843 #ifdef VL_COMPLEX
00844
00845 #include "vl/VLc.h"
00846 Void ComplexTest()
00847 {
00848 cout << "\n+ Test Complex\n\n";
00849
00850 SparseMatc a(4, 6);
00851 SparseVecc b(6);
00852
00853 a.MakeDiag(float_complex(0.5, 0.8) * 2.5);
00854 cout << "a:\n" << a << endl;
00855 b = vl_y;
00856 cout << "b: " << b << endl;
00857
00858 cout << "a * b = " << a * b << endl;
00859
00860 SparseMatc na(4, 6);
00861 na.MakeDiag(5.1);
00862
00863 cout << "na:\n" << na << endl;
00864 cout << "na * b = " << na * b << endl;
00865
00866 a = vl_Z;
00867 b = vl_0;
00868 a = vl_I;
00869 b = vl_x;
00870
00871 Mat4c c;
00872 Vec4c d;
00873
00874 c = vl_1;
00875 d = vl_y;
00876
00877 cout << "c:\n" << c << endl;
00878 cout << "d:" << d << endl;
00879
00880 cout << "ca * cb = " << c * d + d << endl;
00881 }
00882 #endif
00883
