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


kni
Author(s): Martin Günther
autogenerated on Mon Aug 14 2017 02:44:13