tmte.cc
Go to the documentation of this file.
00001 
00002 //#define WANT_STREAM
00003 #define WANT_MATH
00004 
00005 #include "include.h"
00006 
00007 #include "newmatap.h"
00008 //#include "newmatio.h"
00009 
00010 #include "tmt.h"
00011 
00012 #ifdef use_namespace
00013 using namespace NEWMAT;
00014 #endif
00015 
00016 // check D is sorted
00017 void CheckIsSorted(const DiagonalMatrix& D, bool ascending = false)
00018 {
00019    DiagonalMatrix D1 = D;
00020    if (ascending) SortAscending(D1); else SortDescending(D1);
00021    D1 -= D; Print(D1);
00022 }
00023 
00024 
00025 
00026 void trymate()
00027 {
00028 
00029 
00030    Tracer et("Fourteenth test of Matrix package");
00031    Tracer::PrintTrace();
00032 
00033    {
00034       Tracer et1("Stage 1");
00035       Matrix A(8,5);
00036       {
00037 #ifndef ATandT
00038          Real   a[] = { 22, 10,  2,  3,  7,
00039                         14,  7, 10,  0,  8,
00040                         -1, 13, -1,-11,  3,
00041                         -3, -2, 13, -2,  4,
00042                          9,  8,  1, -2,  4,
00043                          9,  1, -7,  5, -1,
00044                          2, -6,  6,  5,  1,
00045                          4,  5,  0, -2,  2 };
00046 #else
00047          Real a[40];
00048          a[ 0]=22; a[ 1]=10; a[ 2]= 2; a[ 3]= 3; a[ 4]= 7;
00049          a[ 5]=14; a[ 6]= 7; a[ 7]=10; a[ 8]= 0; a[ 9]= 8;
00050          a[10]=-1; a[11]=13; a[12]=-1; a[13]=-11;a[14]= 3;
00051          a[15]=-3; a[16]=-2; a[17]=13; a[18]=-2; a[19]= 4;
00052          a[20]= 9; a[21]= 8; a[22]= 1; a[23]=-2; a[24]= 4;
00053          a[25]= 9; a[26]= 1; a[27]=-7; a[28]= 5; a[29]=-1;
00054          a[30]= 2; a[31]=-6; a[32]= 6; a[33]= 5; a[34]= 1;
00055          a[35]= 4; a[36]= 5; a[37]= 0; a[38]=-2; a[39]= 2;
00056 #endif
00057          A << a;
00058       }
00059       DiagonalMatrix D; Matrix U; Matrix V;
00060       int anc = A.Ncols(); IdentityMatrix I(anc);
00061       SymmetricMatrix S1; S1 << A.t() * A;
00062       SymmetricMatrix S2; S2 << A * A.t();
00063       Real zero = 0.0; SVD(A+zero,D,U,V); CheckIsSorted(D);
00064       DiagonalMatrix D1; SVD(A,D1); CheckIsSorted(D1);
00065       D1 -= D; Clean(D1,0.000000001);Print(D1);
00066       Matrix W;
00067       SVD(A, D1, W, W, true, false); D1 -= D; W -= U;
00068       Clean(W,0.000000001); Print(W); Clean(D1,0.000000001); Print(D1);
00069       Matrix WX;
00070       SVD(A, D1, WX, W, false, true); D1 -= D; W -= V;
00071       Clean(W,0.000000001); Print(W); Clean(D1,0.000000001); Print(D1);
00072       Matrix SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
00073       Matrix SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
00074       Matrix B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
00075       D1=0.0;  SVD(A,D1,A); CheckIsSorted(D1);
00076       A -= U; Clean(A,0.000000001); Print(A);
00077       D(1) -= sqrt(1248.0); D(2) -= 20; D(3) -= sqrt(384.0);
00078       Clean(D,0.000000001); Print(D);
00079 
00080       Jacobi(S1, D, V);  CheckIsSorted(D, true);
00081       V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
00082       D = D.Reverse(); D(1)-=1248; D(2)-=400; D(3)-=384;
00083       Clean(D,0.000000001); Print(D);
00084 
00085       Jacobi(S1, D);  CheckIsSorted(D, true);
00086       D = D.Reverse(); D(1)-=1248; D(2)-=400; D(3)-=384;
00087       Clean(D,0.000000001); Print(D);
00088 
00089       SymmetricMatrix JW(5);
00090       Jacobi(S1, D, JW);  CheckIsSorted(D, true);
00091       D = D.Reverse(); D(1)-=1248; D(2)-=400; D(3)-=384;
00092       Clean(D,0.000000001); Print(D);
00093 
00094       Jacobi(S2, D, V);  CheckIsSorted(D, true);
00095       V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
00096       D = D.Reverse(); D(1)-=1248; D(2)-=400; D(3)-=384;
00097       Clean(D,0.000000001); Print(D);
00098 
00099       EigenValues(S1, D, V); CheckIsSorted(D, true);
00100       V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
00101       D(5)-=1248; D(4)-=400; D(3)-=384;
00102       Clean(D,0.000000001); Print(D);
00103 
00104       EigenValues(S2, D, V); CheckIsSorted(D, true);
00105       V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
00106       D(8)-=1248; D(7)-=400; D(6)-=384;
00107       Clean(D,0.000000001); Print(D);
00108 
00109       EigenValues(S1, D); CheckIsSorted(D, true);
00110       D(5)-=1248; D(4)-=400; D(3)-=384;
00111       Clean(D,0.000000001); Print(D);
00112 
00113       SymmetricMatrix EW(S2);
00114       EigenValues(S2, D, EW); CheckIsSorted(D, true);
00115       D(8)-=1248; D(7)-=400; D(6)-=384;
00116       Clean(D,0.000000001); Print(D);
00117 
00118    }
00119 
00120    {
00121       Tracer et1("Stage 2");
00122       Matrix A(20,21);
00123       int i,j;
00124       for (i=1; i<=20; i++) for (j=1; j<=21; j++)
00125       { if (i>j) A(i,j) = 0; else if (i==j) A(i,j) = 21-i; else A(i,j) = -1; }
00126       A = A.t();
00127       SymmetricMatrix S1; S1 << A.t() * A;
00128       SymmetricMatrix S2; S2 << A * A.t();
00129       DiagonalMatrix D; Matrix U; Matrix V;
00130 #ifdef ATandT
00131       int anc = A.Ncols(); DiagonalMatrix I(anc);     // AT&T 2.1 bug
00132 #else
00133       DiagonalMatrix I(A.Ncols());
00134 #endif
00135       I=1.0;
00136       SVD(A,D,U,V); CheckIsSorted(D);
00137       Matrix SU = U.t() * U - I;    Clean(SU,0.000000001); Print(SU);
00138       Matrix SV = V.t() * V - I;    Clean(SV,0.000000001); Print(SV);
00139       Matrix B = U * D * V.t() - A; Clean(B,0.000000001);  Print(B);
00140       for (i=1; i<=20; i++)  D(i) -= sqrt((22.0-i)*(21.0-i));
00141       Clean(D,0.000000001); Print(D);
00142       Jacobi(S1, D, V); CheckIsSorted(D, true);
00143       V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
00144       D = D.Reverse();
00145       for (i=1; i<=20; i++)  D(i) -= (22-i)*(21-i);
00146       Clean(D,0.000000001); Print(D);
00147       Jacobi(S2, D, V); CheckIsSorted(D, true);
00148       V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
00149       D = D.Reverse();
00150       for (i=1; i<=20; i++)  D(i) -= (22-i)*(21-i);
00151       Clean(D,0.000000001); Print(D);
00152 
00153       EigenValues(S1, D, V); CheckIsSorted(D, true);
00154       V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
00155       for (i=1; i<=20; i++)  D(i) -= (i+1)*i;
00156       Clean(D,0.000000001); Print(D);
00157       EigenValues(S2, D, V); CheckIsSorted(D, true);
00158       V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
00159       for (i=2; i<=21; i++)  D(i) -= (i-1)*i;
00160       Clean(D,0.000000001); Print(D);
00161 
00162       EigenValues(S1, D); CheckIsSorted(D, true);
00163       for (i=1; i<=20; i++)  D(i) -= (i+1)*i;
00164       Clean(D,0.000000001); Print(D);
00165       EigenValues(S2, D); CheckIsSorted(D, true);
00166       for (i=2; i<=21; i++)  D(i) -= (i-1)*i;
00167       Clean(D,0.000000001); Print(D);
00168    }
00169 
00170    {
00171       Tracer et1("Stage 3");
00172       Matrix A(30,30);
00173       int i,j;
00174       for (i=1; i<=30; i++) for (j=1; j<=30; j++)
00175       { if (i>j) A(i,j) = 0; else if (i==j) A(i,j) = 1; else A(i,j) = -1; }
00176       Real d1 = A.LogDeterminant().Value();
00177       DiagonalMatrix D; Matrix U; Matrix V;
00178 #ifdef ATandT
00179       int anc = A.Ncols(); DiagonalMatrix I(anc);     // AT&T 2.1 bug
00180 #else
00181       DiagonalMatrix I(A.Ncols());
00182 #endif
00183       I=1.0;
00184       SVD(A,D,U,V); CheckIsSorted(D);
00185       Matrix SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
00186       Matrix SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
00187       Real d2 = D.LogDeterminant().Value();
00188       Matrix B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
00189       Real d3 = D.LogDeterminant().Value();
00190       ColumnVector Test(3);
00191       Test(1) = d1 - 1; Test(2) = d2 - 1; Test(3) = d3 - 1;
00192       Clean(Test,0.00000001); Print(Test); // only 8 decimal figures
00193       A.ReSize(2,2);
00194       Real a = 1.5; Real b = 2; Real c = 2 * (a*a + b*b);
00195       A << a << b << a << b;
00196       I.ReSize(2); I=1;
00197       SVD(A,D,U,V); CheckIsSorted(D);
00198       SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
00199       SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
00200       B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
00201       D = D*D; SortDescending(D);
00202       DiagonalMatrix D50(2); D50 << c << 0; D = D - D50;
00203       Clean(D,0.000000001);
00204       Print(D);
00205       A << a << a << b << b;
00206       SVD(A,D,U,V); CheckIsSorted(D);
00207       SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
00208       SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
00209       B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
00210       D = D*D; SortDescending(D);
00211       D = D - D50;
00212       Clean(D,0.000000001);
00213       Print(D);
00214    }
00215 
00216    {
00217       Tracer et1("Stage 4");
00218 
00219       // test for bug found by Olof Runborg,
00220       // Department of Numerical Analysis and Computer Science (NADA),
00221       // KTH, Stockholm
00222 
00223       Matrix A(22,20);
00224 
00225       A = 0;
00226 
00227       int a=1;
00228 
00229       A(a+0,a+2) = 1;     A(a+0,a+18) = -1;
00230       A(a+1,a+9) = 1;     A(a+1,a+12) = -1;
00231       A(a+2,a+11) = 1;    A(a+2,a+12) = -1;
00232       A(a+3,a+10) = 1;    A(a+3,a+19) = -1;
00233       A(a+4,a+16) = 1;    A(a+4,a+19) = -1;
00234       A(a+5,a+17) = 1;    A(a+5,a+18) = -1;
00235       A(a+6,a+10) = 1;    A(a+6,a+4) = -1;
00236       A(a+7,a+3) = 1;     A(a+7,a+2) = -1;
00237       A(a+8,a+14) = 1;    A(a+8,a+15) = -1;
00238       A(a+9,a+13) = 1;    A(a+9,a+16) = -1;
00239       A(a+10,a+8) = 1;    A(a+10,a+9) = -1;
00240       A(a+11,a+1) = 1;    A(a+11,a+15) = -1;
00241       A(a+12,a+16) = 1;   A(a+12,a+4) = -1;
00242       A(a+13,a+6) = 1;    A(a+13,a+9) = -1;
00243       A(a+14,a+5) = 1;    A(a+14,a+4) = -1;
00244       A(a+15,a+0) = 1;    A(a+15,a+1) = -1;
00245       A(a+16,a+14) = 1;   A(a+16,a+0) = -1;
00246       A(a+17,a+7) = 1;    A(a+17,a+6) = -1;
00247       A(a+18,a+13) = 1;   A(a+18,a+5) = -1;
00248       A(a+19,a+7) = 1;    A(a+19,a+8) = -1;
00249       A(a+20,a+17) = 1;   A(a+20,a+3) = -1;
00250       A(a+21,a+6) = 1;    A(a+21,a+11) = -1;
00251 
00252 
00253       Matrix U, V; DiagonalMatrix S;
00254 
00255       SVD(A, S, U, V, true, true); CheckIsSorted(S);
00256 
00257       DiagonalMatrix D(20); D = 1;
00258 
00259       Matrix tmp = U.t() * U - D;
00260       Clean(tmp,0.000000001); Print(tmp);
00261 
00262       tmp = V.t() * V - D;
00263       Clean(tmp,0.000000001); Print(tmp);
00264 
00265       tmp = U * S * V.t() - A ;
00266       Clean(tmp,0.000000001); Print(tmp);
00267 
00268    }
00269 
00270    {
00271       Tracer et1("Stage 5");
00272       Matrix A(10,10);
00273 
00274       A.Row(1)  <<  1.00 <<  0.07 <<  0.05 <<  0.00 <<  0.06
00275                 <<  0.09 <<  0.03 <<  0.02 <<  0.02 << -0.03;
00276       A.Row(2)  <<  0.07 <<  1.00 <<  0.05 <<  0.05 << -0.03
00277                 <<  0.07 <<  0.00 <<  0.07 <<  0.00 <<  0.02;
00278       A.Row(3)  <<  0.05 <<  0.05 <<  1.00 <<  0.05 <<  0.02
00279                 <<  0.01 << -0.05 <<  0.04 <<  0.05 << -0.03;
00280       A.Row(4)  <<  0.00 <<  0.05 <<  0.05 <<  1.00 << -0.05
00281                 <<  0.04 <<  0.01 <<  0.02 << -0.05 <<  0.00;
00282       A.Row(5)  <<  0.06 << -0.03 <<  0.02 << -0.05 <<  1.00
00283                 << -0.03 <<  0.02 << -0.02 <<  0.04 <<  0.00;
00284       A.Row(6)  <<  0.09 <<  0.07 <<  0.01 <<  0.04 << -0.03
00285                 <<  1.00 << -0.06 <<  0.08 << -0.02 << -0.10;
00286       A.Row(7)  <<  0.03 <<  0.00 << -0.05 <<  0.01 <<  0.02
00287                 << -0.06 <<  1.00 <<  0.09 <<  0.12 << -0.03;
00288       A.Row(8)  <<  0.02 <<  0.07 <<  0.04 <<  0.02 << -0.02
00289                 <<  0.08 <<  0.09 <<  1.00 <<  0.00 << -0.02;
00290       A.Row(9)  <<  0.02 <<  0.00 <<  0.05 << -0.05 <<  0.04
00291                 << -0.02 <<  0.12 <<  0.00 <<  1.00 <<  0.02;
00292       A.Row(10) << -0.03 <<  0.02 << -0.03 <<  0.00 <<  0.00
00293                 << -0.10 << -0.03 << -0.02 <<  0.02 <<  1.00;
00294 
00295       SymmetricMatrix AS; AS << A;
00296       Matrix V; DiagonalMatrix D, D1;
00297       ColumnVector Check(6);
00298       EigenValues(AS,D,V); CheckIsSorted(D, true);
00299       Check(1) = MaximumAbsoluteValue(A - V * D * V.t());
00300       DiagonalMatrix I(10); I = 1;
00301       Check(2) = MaximumAbsoluteValue(V * V.t() - I);
00302       Check(3) = MaximumAbsoluteValue(V.t() * V - I);
00303 
00304       EigenValues(AS, D1); CheckIsSorted(D1, true);
00305       D -= D1;
00306       Clean(D,0.000000001); Print(D);
00307 
00308       Jacobi(AS,D,V);
00309       Check(4) = MaximumAbsoluteValue(A - V * D * V.t());
00310       Check(5) = MaximumAbsoluteValue(V * V.t() - I);
00311       Check(6) = MaximumAbsoluteValue(V.t() * V - I);
00312 
00313       SortAscending(D);
00314       D -= D1;
00315       Clean(D,0.000000001); Print(D);
00316 
00317       Clean(Check,0.000000001); Print(Check);
00318 
00319       // Check loading rows
00320 
00321       SymmetricMatrix B(10);
00322 
00323       B.Row(1)  <<  1.00;
00324       B.Row(2)  <<  0.07 <<  1.00;
00325       B.Row(3)  <<  0.05 <<  0.05 <<  1.00;
00326       B.Row(4)  <<  0.00 <<  0.05 <<  0.05 <<  1.00;
00327       B.Row(5)  <<  0.06 << -0.03 <<  0.02 << -0.05 <<  1.00;
00328       B.Row(6)  <<  0.09 <<  0.07 <<  0.01 <<  0.04 << -0.03
00329                 <<  1.00;
00330       B.Row(7)  <<  0.03 <<  0.00 << -0.05 <<  0.01 <<  0.02
00331                 << -0.06 <<  1.00;
00332       B.Row(8)  <<  0.02 <<  0.07 <<  0.04 <<  0.02 << -0.02
00333                 <<  0.08 <<  0.09 <<  1.00;
00334       B.Row(9)  <<  0.02 <<  0.00 <<  0.05 << -0.05 <<  0.04
00335                 << -0.02 <<  0.12 <<  0.00 <<  1.00;
00336       B.Row(10) << -0.03 <<  0.02 << -0.03 <<  0.00 <<  0.00
00337                 << -0.10 << -0.03 << -0.02 <<  0.02 <<  1.00;
00338 
00339       B -= AS; Print(B);
00340 
00341    }
00342 
00343    {
00344       Tracer et1("Stage 6");
00345       // badly scaled matrix
00346       Matrix A(9,9);
00347 
00348       A.Row(1) << 1.13324e+012 << 3.68788e+011 << 3.35163e+009
00349                << 3.50193e+011 << 1.25335e+011 << 1.02212e+009
00350                << 3.16602e+009 << 1.02418e+009 << 9.42959e+006;
00351       A.Row(2) << 3.68788e+011 << 1.67128e+011 << 1.27449e+009
00352                << 1.25335e+011 << 6.05413e+010 << 4.34573e+008
00353                << 1.02418e+009 << 4.69192e+008 << 3.61098e+006;
00354       A.Row(3) << 3.35163e+009 << 1.27449e+009 << 1.25571e+007
00355                << 1.02212e+009 << 4.34573e+008 << 3.69769e+006
00356                << 9.42959e+006 << 3.61098e+006 << 3.59450e+004;
00357       A.Row(4) << 3.50193e+011 << 1.25335e+011 << 1.02212e+009
00358                << 1.43514e+011 << 5.42310e+010 << 4.15822e+008
00359                << 1.23068e+009 << 4.31545e+008 << 3.58714e+006;
00360       A.Row(5) << 1.25335e+011 << 6.05413e+010 << 4.34573e+008
00361                << 5.42310e+010 << 2.76601e+010 << 1.89102e+008
00362                << 4.31545e+008 << 2.09778e+008 << 1.51083e+006;
00363       A.Row(6) << 1.02212e+009 << 4.34573e+008 << 3.69769e+006
00364                << 4.15822e+008 << 1.89102e+008 << 1.47143e+006
00365                << 3.58714e+006 << 1.51083e+006 << 1.30165e+004;
00366       A.Row(7) << 3.16602e+009 << 1.02418e+009 << 9.42959e+006
00367                << 1.23068e+009 << 4.31545e+008 << 3.58714e+006
00368                << 1.12335e+007 << 3.54778e+006 << 3.34311e+004;
00369       A.Row(8) << 1.02418e+009 << 4.69192e+008 << 3.61098e+006
00370                << 4.31545e+008 << 2.09778e+008 << 1.51083e+006
00371                << 3.54778e+006 << 1.62552e+006 << 1.25885e+004;
00372       A.Row(9) << 9.42959e+006 << 3.61098e+006 << 3.59450e+004
00373                << 3.58714e+006 << 1.51083e+006 << 1.30165e+004
00374                << 3.34311e+004 << 1.25885e+004 << 1.28000e+002;
00375 
00376 
00377       SymmetricMatrix AS; AS << A;
00378       Matrix V; DiagonalMatrix D, D1;
00379       ColumnVector Check(6);
00380       EigenValues(AS,D,V); CheckIsSorted(D, true);
00381       Check(1) = MaximumAbsoluteValue(A - V * D * V.t()) / 100000;
00382       DiagonalMatrix I(9); I = 1;
00383       Check(2) = MaximumAbsoluteValue(V * V.t() - I);
00384       Check(3) = MaximumAbsoluteValue(V.t() * V - I);
00385 
00386       EigenValues(AS, D1);
00387       D -= D1;
00388       Clean(D,0.001); Print(D);
00389 
00390       Jacobi(AS,D,V);
00391       Check(4) = MaximumAbsoluteValue(A - V * D * V.t()) / 100000;
00392       Check(5) = MaximumAbsoluteValue(V * V.t() - I);
00393       Check(6) = MaximumAbsoluteValue(V.t() * V - I);
00394 
00395       SortAscending(D);
00396       D -= D1;
00397       Clean(D,0.001); Print(D);
00398 
00399       Clean(Check,0.0000001); Print(Check);
00400    }
00401 
00402    {
00403       Tracer et1("Stage 7");
00404       // matrix with all singular values close to 1
00405       Matrix A(8,8);
00406       A.Row(1)<<-0.4343<<-0.0445<<-0.4582<<-0.1612<<-0.3191<<-0.6784<<0.1068<<0;
00407       A.Row(2)<<0.5791<<0.5517<<0.2575<<-0.1055<<-0.0437<<-0.5282<<0.0442<<0;
00408       A.Row(3)<<0.5709<<-0.5179<<-0.3275<<0.2598<<-0.196<<-0.1451<<-0.4143<<0;
00409       A.Row(4)<<0.2785<<-0.5258<<0.1251<<-0.4382<<0.0514<<-0.0446<<0.6586<<0;
00410       A.Row(5)<<0.2654<<0.3736<<-0.7436<<-0.0122<<0.0376<<0.3465<<0.3397<<0;
00411       A.Row(6)<<0.0173<<-0.0056<<-0.1903<<-0.7027<<0.4863<<-0.0199<<-0.4825<<0;
00412       A.Row(7)<<0.0434<<0.0966<<0.1083<<-0.4576<<-0.7857<<0.3425<<-0.1818<<0;
00413       A.Row(8)<<0.0<<0.0<<0.0<<0.0<<0.0<<0.0<<0.0<<-1.0;
00414       Matrix U,V; DiagonalMatrix D;
00415       SVD(A,D,U,V); CheckIsSorted(D);
00416       Matrix B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
00417       DiagonalMatrix I(8); I = 1; D -= I; Clean(D,0.0001); Print(D);
00418       U *= U.t(); U -= I; Clean(U,0.000000001); Print(U);
00419       V *= V.t(); V -= I; Clean(V,0.000000001); Print(V);
00420 
00421    }
00422 
00423    {
00424       Tracer et1("Stage 8");
00425       // check SortSV functions
00426 
00427       Matrix A(15, 10);
00428       int i, j;
00429       for (i = 1; i <= 15; ++i) for (j = 1; j <= 10; ++j)
00430          A(i, j) = i + j / 1000.0;
00431       DiagonalMatrix D(10);
00432       D << 0.2 << 0.5 << 0.1 << 0.7 << 0.8 << 0.3 << 0.4 << 0.7 << 0.9 << 0.6;
00433       Matrix U = A; Matrix V = 10 - 2 * A;
00434       Matrix Prod = U * D * V.t();
00435 
00436       DiagonalMatrix D2 = D; SortDescending(D2);
00437       DiagonalMatrix D1 = D; SortSV(D1, U, V); Matrix X = D1 - D2; Print(X);
00438       X = Prod - U * D1 * V.t(); Clean(X,0.000000001); Print(X);
00439       U = A; V = 10 - 2 * A;
00440       D1 = D; SortSV(D1, U); X = D1 - D2; Print(X);
00441       D1 = D; SortSV(D1, V); X = D1 - D2; Print(X);
00442       X = Prod - U * D1 * V.t(); Clean(X,0.000000001); Print(X);
00443 
00444       D2 = D; SortAscending(D2);
00445       U = A; V = 10 - 2 * A;
00446       D1 = D; SortSV(D1, U, V, true); X = D1 - D2; Print(X);
00447       X = Prod - U * D1 * V.t(); Clean(X,0.000000001); Print(X);
00448       U = A; V = 10 - 2 * A;
00449       D1 = D; SortSV(D1, U, true); X = D1 - D2; Print(X);
00450       D1 = D; SortSV(D1, V, true); X = D1 - D2; Print(X);
00451       X = Prod - U * D1 * V.t(); Clean(X,0.000000001); Print(X);
00452    }
00453 
00454    {
00455       Tracer et1("Stage 9");
00456       // Tom William's example
00457       Matrix A(10,10);
00458       Matrix U;
00459       Matrix V;
00460       DiagonalMatrix Sigma;
00461       Real myVals[] =
00462       {
00463          1,    1,    1,    1,    1,    1,    1,    1,    1,    1,
00464          1,    1,    1,    1,    1,    1,    1,    1,    1,    1,
00465          1,    1,    1,    1,    1,    1,    1,    1,    1,    1,
00466          1,    1,    1,    1,    1,    1,    1,    1,    1,    1,
00467          1,    1,    1,    1,    1,    1,    1,    1,    1,    1,
00468          1,    1,    1,    1,    1,    1,    1,    1,    1,    0,
00469          1,    1,    1,    1,    1,    1,    1,    1,    1,    0,
00470          1,    1,    1,    1,    1,    1,    1,    1,    0,    0,
00471          1,    1,    1,    1,    1,    1,    1,    0,    0,    0,
00472          1,    1,    1,    1,    1,    0,    0,    0,    0,    0,
00473       };
00474 
00475       A << myVals;
00476       SVD(A, Sigma, U, V); CheckIsSorted(Sigma);
00477       A -= U * Sigma * V.t();
00478       Clean(A, 0.000000001); Print(A);
00479    }
00480 
00481 
00482 
00483 }


rl_agent
Author(s): Todd Hester
autogenerated on Thu Jun 6 2019 22:00:13