00001
00002
00003 #define WANT_MATH
00004
00005 #include "include.h"
00006
00007 #include "newmatap.h"
00008
00009
00010 #include "tmt.h"
00011
00012 #ifdef use_namespace
00013 using namespace NEWMAT;
00014 #endif
00015
00016
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);
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);
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);
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
00220
00221
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
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
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
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
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
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 }