Planet
navi homePPSaboutscreenshotsdownloaddevelopmentforum

source: orxonox.OLD/orxonox/trunk/src/lib/newmat/tmte.cpp @ 4608

Last change on this file since 4608 was 4565, checked in by patrick, 19 years ago

orxonox/trunk: added the newmat library to the project. needs some translation in directory, temp under util/newmat. is needed by the collision detection engine to perform lin alg operations such as eigenvector decomposition. perhaps we will make our own library to do that later.
File size: 17.4 KB
Line 
1
2//#define WANT_STREAM
3#define WANT_MATH
4
5#include "include.h"
6
7#include "newmatap.h"
8//#include "newmatio.h"
9
10#include "tmt.h"
11
12#ifdef use_namespace
13using namespace NEWMAT;
14#endif
15
16// check D is sorted
17void CheckIsSorted(const DiagonalMatrix& D, bool ascending = false)
18{
19   DiagonalMatrix D1 = D;
20   if (ascending) SortAscending(D1); else SortDescending(D1);
21   D1 -= D; Print(D1);
22}
23
24
25
26void trymate()
27{
28
29
30   Tracer et("Fourteenth test of Matrix package");
31   Tracer::PrintTrace();
32
33   {
34      Tracer et1("Stage 1");
35      Matrix A(8,5);
36      {
37#ifndef ATandT
38         Real   a[] = { 22, 10,  2,  3,  7,
39                        14,  7, 10,  0,  8,
40                        -1, 13, -1,-11,  3,
41                        -3, -2, 13, -2,  4,
42                         9,  8,  1, -2,  4,
43                         9,  1, -7,  5, -1,
44                         2, -6,  6,  5,  1,
45                         4,  5,  0, -2,  2 };
46#else
47         Real a[40];
48         a[ 0]=22; a[ 1]=10; a[ 2]= 2; a[ 3]= 3; a[ 4]= 7;
49         a[ 5]=14; a[ 6]= 7; a[ 7]=10; a[ 8]= 0; a[ 9]= 8;
50         a[10]=-1; a[11]=13; a[12]=-1; a[13]=-11;a[14]= 3;
51         a[15]=-3; a[16]=-2; a[17]=13; a[18]=-2; a[19]= 4;
52         a[20]= 9; a[21]= 8; a[22]= 1; a[23]=-2; a[24]= 4;
53         a[25]= 9; a[26]= 1; a[27]=-7; a[28]= 5; a[29]=-1;
54         a[30]= 2; a[31]=-6; a[32]= 6; a[33]= 5; a[34]= 1;
55         a[35]= 4; a[36]= 5; a[37]= 0; a[38]=-2; a[39]= 2;
56#endif
57         A << a;
58      }
59      DiagonalMatrix D; Matrix U; Matrix V;
60      int anc = A.Ncols(); IdentityMatrix I(anc);
61      SymmetricMatrix S1; S1 << A.t() * A;
62      SymmetricMatrix S2; S2 << A * A.t();
63      Real zero = 0.0; SVD(A+zero,D,U,V); CheckIsSorted(D);
64      DiagonalMatrix D1; SVD(A,D1); CheckIsSorted(D1);
65      Print(DiagonalMatrix(D-D1));
66      Matrix W;
67      SVD(A, D1, W, W, true, false); D1 -= D; W -= U;
68      Clean(W,0.000000001); Print(W); Clean(D1,0.000000001); Print(D1);
69      Matrix WX;
70      SVD(A, D1, WX, W, false, true); D1 -= D; W -= V;
71      Clean(W,0.000000001); Print(W); Clean(D1,0.000000001); Print(D1);
72      Matrix SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
73      Matrix SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
74      Matrix B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
75      D1=0.0;  SVD(A,D1,A); CheckIsSorted(D1); Print(Matrix(A-U));
76      D(1) -= sqrt(1248.0); D(2) -= 20; D(3) -= sqrt(384.0);
77      Clean(D,0.000000001); Print(D);
78
79      Jacobi(S1, D, V);  CheckIsSorted(D, true);
80      V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
81      D = D.Reverse(); D(1)-=1248; D(2)-=400; D(3)-=384;
82      Clean(D,0.000000001); Print(D);
83
84      Jacobi(S1, D);  CheckIsSorted(D, true);
85      D = D.Reverse(); D(1)-=1248; D(2)-=400; D(3)-=384;
86      Clean(D,0.000000001); Print(D);
87
88      SymmetricMatrix JW(5);
89      Jacobi(S1, D, JW);  CheckIsSorted(D, true);
90      D = D.Reverse(); D(1)-=1248; D(2)-=400; D(3)-=384;
91      Clean(D,0.000000001); Print(D);
92
93      Jacobi(S2, D, V);  CheckIsSorted(D, true);
94      V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
95      D = D.Reverse(); D(1)-=1248; D(2)-=400; D(3)-=384;
96      Clean(D,0.000000001); Print(D);
97
98      EigenValues(S1, D, V); CheckIsSorted(D, true);
99      V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
100      D(5)-=1248; D(4)-=400; D(3)-=384;
101      Clean(D,0.000000001); Print(D);
102
103      EigenValues(S2, D, V); CheckIsSorted(D, true);
104      V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
105      D(8)-=1248; D(7)-=400; D(6)-=384;
106      Clean(D,0.000000001); Print(D);
107
108      EigenValues(S1, D); CheckIsSorted(D, true);
109      D(5)-=1248; D(4)-=400; D(3)-=384;
110      Clean(D,0.000000001); Print(D);
111
112      SymmetricMatrix EW(S2);
113      EigenValues(S2, D, EW); CheckIsSorted(D, true);
114      D(8)-=1248; D(7)-=400; D(6)-=384;
115      Clean(D,0.000000001); Print(D);
116
117   }
118
119   {
120      Tracer et1("Stage 2");
121      Matrix A(20,21);
122      int i,j;
123      for (i=1; i<=20; i++) for (j=1; j<=21; j++)
124      { if (i>j) A(i,j) = 0; else if (i==j) A(i,j) = 21-i; else A(i,j) = -1; }
125      A = A.t();
126      SymmetricMatrix S1; S1 << A.t() * A;
127      SymmetricMatrix S2; S2 << A * A.t();
128      DiagonalMatrix D; Matrix U; Matrix V;
129#ifdef ATandT
130      int anc = A.Ncols(); DiagonalMatrix I(anc);     // AT&T 2.1 bug
131#else
132      DiagonalMatrix I(A.Ncols());
133#endif
134      I=1.0;
135      SVD(A,D,U,V); CheckIsSorted(D);
136      Matrix SU = U.t() * U - I;    Clean(SU,0.000000001); Print(SU);
137      Matrix SV = V.t() * V - I;    Clean(SV,0.000000001); Print(SV);
138      Matrix B = U * D * V.t() - A; Clean(B,0.000000001);  Print(B);
139      for (i=1; i<=20; i++)  D(i) -= sqrt((22.0-i)*(21.0-i));
140      Clean(D,0.000000001); Print(D);
141      Jacobi(S1, D, V); CheckIsSorted(D, true);
142      V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
143      D = D.Reverse();
144      for (i=1; i<=20; i++)  D(i) -= (22-i)*(21-i);
145      Clean(D,0.000000001); Print(D);
146      Jacobi(S2, D, V); CheckIsSorted(D, true);
147      V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
148      D = D.Reverse();
149      for (i=1; i<=20; i++)  D(i) -= (22-i)*(21-i);
150      Clean(D,0.000000001); Print(D);
151
152      EigenValues(S1, D, V); CheckIsSorted(D, true);
153      V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
154      for (i=1; i<=20; i++)  D(i) -= (i+1)*i;
155      Clean(D,0.000000001); Print(D);
156      EigenValues(S2, D, V); CheckIsSorted(D, true);
157      V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
158      for (i=2; i<=21; i++)  D(i) -= (i-1)*i;
159      Clean(D,0.000000001); Print(D);
160
161      EigenValues(S1, D); CheckIsSorted(D, true);
162      for (i=1; i<=20; i++)  D(i) -= (i+1)*i;
163      Clean(D,0.000000001); Print(D);
164      EigenValues(S2, D); CheckIsSorted(D, true);
165      for (i=2; i<=21; i++)  D(i) -= (i-1)*i;
166      Clean(D,0.000000001); Print(D);
167   }
168
169   {
170      Tracer et1("Stage 3");
171      Matrix A(30,30);
172      int i,j;
173      for (i=1; i<=30; i++) for (j=1; j<=30; j++)
174      { if (i>j) A(i,j) = 0; else if (i==j) A(i,j) = 1; else A(i,j) = -1; }
175      Real d1 = A.LogDeterminant().Value();
176      DiagonalMatrix D; Matrix U; Matrix V;
177#ifdef ATandT
178      int anc = A.Ncols(); DiagonalMatrix I(anc);     // AT&T 2.1 bug
179#else
180      DiagonalMatrix I(A.Ncols());
181#endif
182      I=1.0;
183      SVD(A,D,U,V); CheckIsSorted(D);
184      Matrix SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
185      Matrix SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
186      Real d2 = D.LogDeterminant().Value();
187      Matrix B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
188      Real d3 = D.LogDeterminant().Value();
189      ColumnVector Test(3);
190      Test(1) = d1 - 1; Test(2) = d2 - 1; Test(3) = d3 - 1;
191      Clean(Test,0.00000001); Print(Test); // only 8 decimal figures
192      A.ReSize(2,2);
193      Real a = 1.5; Real b = 2; Real c = 2 * (a*a + b*b);
194      A << a << b << a << b;
195      I.ReSize(2); I=1;
196      SVD(A,D,U,V); CheckIsSorted(D);
197      SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
198      SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
199      B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
200      D = D*D; SortDescending(D);
201      DiagonalMatrix D50(2); D50 << c << 0; D = D - D50;
202      Clean(D,0.000000001);
203      Print(D);
204      A << a << a << b << b;
205      SVD(A,D,U,V); CheckIsSorted(D);
206      SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
207      SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
208      B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
209      D = D*D; SortDescending(D);
210      D = D - D50;
211      Clean(D,0.000000001);
212      Print(D);
213   }
214
215   {
216      Tracer et1("Stage 4");
217
218      // test for bug found by Olof Runborg,
219      // Department of Numerical Analysis and Computer Science (NADA),
220      // KTH, Stockholm
221
222      Matrix A(22,20);
223
224      A = 0;
225
226      int a=1;
227
228      A(a+0,a+2) = 1;     A(a+0,a+18) = -1;
229      A(a+1,a+9) = 1;     A(a+1,a+12) = -1;
230      A(a+2,a+11) = 1;    A(a+2,a+12) = -1;
231      A(a+3,a+10) = 1;    A(a+3,a+19) = -1;
232      A(a+4,a+16) = 1;    A(a+4,a+19) = -1;
233      A(a+5,a+17) = 1;    A(a+5,a+18) = -1;
234      A(a+6,a+10) = 1;    A(a+6,a+4) = -1;
235      A(a+7,a+3) = 1;     A(a+7,a+2) = -1;
236      A(a+8,a+14) = 1;    A(a+8,a+15) = -1;
237      A(a+9,a+13) = 1;    A(a+9,a+16) = -1;
238      A(a+10,a+8) = 1;    A(a+10,a+9) = -1;
239      A(a+11,a+1) = 1;    A(a+11,a+15) = -1;
240      A(a+12,a+16) = 1;   A(a+12,a+4) = -1;
241      A(a+13,a+6) = 1;    A(a+13,a+9) = -1;
242      A(a+14,a+5) = 1;    A(a+14,a+4) = -1;
243      A(a+15,a+0) = 1;    A(a+15,a+1) = -1;
244      A(a+16,a+14) = 1;   A(a+16,a+0) = -1;
245      A(a+17,a+7) = 1;    A(a+17,a+6) = -1;
246      A(a+18,a+13) = 1;   A(a+18,a+5) = -1;
247      A(a+19,a+7) = 1;    A(a+19,a+8) = -1;
248      A(a+20,a+17) = 1;   A(a+20,a+3) = -1;
249      A(a+21,a+6) = 1;    A(a+21,a+11) = -1;
250
251
252      Matrix U, V; DiagonalMatrix S;
253
254      SVD(A, S, U, V, true, true); CheckIsSorted(S);
255
256      DiagonalMatrix D(20); D = 1;
257
258      Matrix tmp = U.t() * U - D;
259      Clean(tmp,0.000000001); Print(tmp);
260
261      tmp = V.t() * V - D;
262      Clean(tmp,0.000000001); Print(tmp);
263
264      tmp = U * S * V.t() - A ;
265      Clean(tmp,0.000000001); Print(tmp);
266
267   }
268
269   {
270      Tracer et1("Stage 5");
271      Matrix A(10,10);
272
273      A.Row(1)  <<  1.00 <<  0.07 <<  0.05 <<  0.00 <<  0.06
274                <<  0.09 <<  0.03 <<  0.02 <<  0.02 << -0.03;
275      A.Row(2)  <<  0.07 <<  1.00 <<  0.05 <<  0.05 << -0.03
276                <<  0.07 <<  0.00 <<  0.07 <<  0.00 <<  0.02;
277      A.Row(3)  <<  0.05 <<  0.05 <<  1.00 <<  0.05 <<  0.02
278                <<  0.01 << -0.05 <<  0.04 <<  0.05 << -0.03;
279      A.Row(4)  <<  0.00 <<  0.05 <<  0.05 <<  1.00 << -0.05
280                <<  0.04 <<  0.01 <<  0.02 << -0.05 <<  0.00;
281      A.Row(5)  <<  0.06 << -0.03 <<  0.02 << -0.05 <<  1.00
282                << -0.03 <<  0.02 << -0.02 <<  0.04 <<  0.00;
283      A.Row(6)  <<  0.09 <<  0.07 <<  0.01 <<  0.04 << -0.03
284                <<  1.00 << -0.06 <<  0.08 << -0.02 << -0.10;
285      A.Row(7)  <<  0.03 <<  0.00 << -0.05 <<  0.01 <<  0.02
286                << -0.06 <<  1.00 <<  0.09 <<  0.12 << -0.03;
287      A.Row(8)  <<  0.02 <<  0.07 <<  0.04 <<  0.02 << -0.02
288                <<  0.08 <<  0.09 <<  1.00 <<  0.00 << -0.02;
289      A.Row(9)  <<  0.02 <<  0.00 <<  0.05 << -0.05 <<  0.04
290                << -0.02 <<  0.12 <<  0.00 <<  1.00 <<  0.02;
291      A.Row(10) << -0.03 <<  0.02 << -0.03 <<  0.00 <<  0.00
292                << -0.10 << -0.03 << -0.02 <<  0.02 <<  1.00;
293
294      SymmetricMatrix AS; AS << A;
295      Matrix V; DiagonalMatrix D, D1;
296      ColumnVector Check(6);
297      EigenValues(AS,D,V); CheckIsSorted(D, true);
298      Check(1) = MaximumAbsoluteValue(A - V * D * V.t());
299      DiagonalMatrix I(10); I = 1;
300      Check(2) = MaximumAbsoluteValue(V * V.t() - I);
301      Check(3) = MaximumAbsoluteValue(V.t() * V - I);
302
303      EigenValues(AS, D1); CheckIsSorted(D1, true);
304      D -= D1;
305      Clean(D,0.000000001); Print(D);
306
307      Jacobi(AS,D,V);
308      Check(4) = MaximumAbsoluteValue(A - V * D * V.t());
309      Check(5) = MaximumAbsoluteValue(V * V.t() - I);
310      Check(6) = MaximumAbsoluteValue(V.t() * V - I);
311
312      SortAscending(D);
313      D -= D1;
314      Clean(D,0.000000001); Print(D);
315
316      Clean(Check,0.000000001); Print(Check);
317
318      // Check loading rows
319
320      SymmetricMatrix B(10);
321
322      B.Row(1)  <<  1.00;
323      B.Row(2)  <<  0.07 <<  1.00;
324      B.Row(3)  <<  0.05 <<  0.05 <<  1.00;
325      B.Row(4)  <<  0.00 <<  0.05 <<  0.05 <<  1.00;
326      B.Row(5)  <<  0.06 << -0.03 <<  0.02 << -0.05 <<  1.00;
327      B.Row(6)  <<  0.09 <<  0.07 <<  0.01 <<  0.04 << -0.03
328                <<  1.00;
329      B.Row(7)  <<  0.03 <<  0.00 << -0.05 <<  0.01 <<  0.02
330                << -0.06 <<  1.00;
331      B.Row(8)  <<  0.02 <<  0.07 <<  0.04 <<  0.02 << -0.02
332                <<  0.08 <<  0.09 <<  1.00;
333      B.Row(9)  <<  0.02 <<  0.00 <<  0.05 << -0.05 <<  0.04
334                << -0.02 <<  0.12 <<  0.00 <<  1.00;
335      B.Row(10) << -0.03 <<  0.02 << -0.03 <<  0.00 <<  0.00
336                << -0.10 << -0.03 << -0.02 <<  0.02 <<  1.00;
337
338      B -= AS; Print(B);
339
340   }
341
342   {
343      Tracer et1("Stage 6");
344      // badly scaled matrix
345      Matrix A(9,9);
346
347      A.Row(1) << 1.13324e+012 << 3.68788e+011 << 3.35163e+009
348               << 3.50193e+011 << 1.25335e+011 << 1.02212e+009
349               << 3.16602e+009 << 1.02418e+009 << 9.42959e+006;
350      A.Row(2) << 3.68788e+011 << 1.67128e+011 << 1.27449e+009
351               << 1.25335e+011 << 6.05413e+010 << 4.34573e+008
352               << 1.02418e+009 << 4.69192e+008 << 3.61098e+006;
353      A.Row(3) << 3.35163e+009 << 1.27449e+009 << 1.25571e+007
354               << 1.02212e+009 << 4.34573e+008 << 3.69769e+006
355               << 9.42959e+006 << 3.61098e+006 << 3.59450e+004;
356      A.Row(4) << 3.50193e+011 << 1.25335e+011 << 1.02212e+009
357               << 1.43514e+011 << 5.42310e+010 << 4.15822e+008
358               << 1.23068e+009 << 4.31545e+008 << 3.58714e+006;
359      A.Row(5) << 1.25335e+011 << 6.05413e+010 << 4.34573e+008
360               << 5.42310e+010 << 2.76601e+010 << 1.89102e+008
361               << 4.31545e+008 << 2.09778e+008 << 1.51083e+006;
362      A.Row(6) << 1.02212e+009 << 4.34573e+008 << 3.69769e+006
363               << 4.15822e+008 << 1.89102e+008 << 1.47143e+006
364               << 3.58714e+006 << 1.51083e+006 << 1.30165e+004;
365      A.Row(7) << 3.16602e+009 << 1.02418e+009 << 9.42959e+006
366               << 1.23068e+009 << 4.31545e+008 << 3.58714e+006
367               << 1.12335e+007 << 3.54778e+006 << 3.34311e+004;
368      A.Row(8) << 1.02418e+009 << 4.69192e+008 << 3.61098e+006
369               << 4.31545e+008 << 2.09778e+008 << 1.51083e+006
370               << 3.54778e+006 << 1.62552e+006 << 1.25885e+004;
371      A.Row(9) << 9.42959e+006 << 3.61098e+006 << 3.59450e+004
372               << 3.58714e+006 << 1.51083e+006 << 1.30165e+004
373               << 3.34311e+004 << 1.25885e+004 << 1.28000e+002;
374
375
376      SymmetricMatrix AS; AS << A;
377      Matrix V; DiagonalMatrix D, D1;
378      ColumnVector Check(6);
379      EigenValues(AS,D,V); CheckIsSorted(D, true);
380      Check(1) = MaximumAbsoluteValue(A - V * D * V.t()) / 100000;
381      DiagonalMatrix I(9); I = 1;
382      Check(2) = MaximumAbsoluteValue(V * V.t() - I);
383      Check(3) = MaximumAbsoluteValue(V.t() * V - I);
384
385      EigenValues(AS, D1);
386      D -= D1;
387      Clean(D,0.001); Print(D);
388
389      Jacobi(AS,D,V);
390      Check(4) = MaximumAbsoluteValue(A - V * D * V.t()) / 100000;
391      Check(5) = MaximumAbsoluteValue(V * V.t() - I);
392      Check(6) = MaximumAbsoluteValue(V.t() * V - I);
393
394      SortAscending(D);
395      D -= D1;
396      Clean(D,0.001); Print(D);
397
398      Clean(Check,0.0000001); Print(Check);
399   }
400
401   {
402      Tracer et1("Stage 7");
403      // matrix with all singular values close to 1
404      Matrix A(8,8);
405      A.Row(1)<<-0.4343<<-0.0445<<-0.4582<<-0.1612<<-0.3191<<-0.6784<<0.1068<<0;
406      A.Row(2)<<0.5791<<0.5517<<0.2575<<-0.1055<<-0.0437<<-0.5282<<0.0442<<0;
407      A.Row(3)<<0.5709<<-0.5179<<-0.3275<<0.2598<<-0.196<<-0.1451<<-0.4143<<0;
408      A.Row(4)<<0.2785<<-0.5258<<0.1251<<-0.4382<<0.0514<<-0.0446<<0.6586<<0;
409      A.Row(5)<<0.2654<<0.3736<<-0.7436<<-0.0122<<0.0376<<0.3465<<0.3397<<0;
410      A.Row(6)<<0.0173<<-0.0056<<-0.1903<<-0.7027<<0.4863<<-0.0199<<-0.4825<<0;
411      A.Row(7)<<0.0434<<0.0966<<0.1083<<-0.4576<<-0.7857<<0.3425<<-0.1818<<0;
412      A.Row(8)<<0.0<<0.0<<0.0<<0.0<<0.0<<0.0<<0.0<<-1.0;
413      Matrix U,V; DiagonalMatrix D;
414      SVD(A,D,U,V); CheckIsSorted(D);
415      Matrix B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
416      DiagonalMatrix I(8); I = 1; D -= I; Clean(D,0.0001); Print(D);
417      U *= U.t(); U -= I; Clean(U,0.000000001); Print(U);
418      V *= V.t(); V -= I; Clean(V,0.000000001); Print(V);
419
420   }
421
422   {
423      Tracer et1("Stage 8");
424      // check SortSV functions
425
426      Matrix A(15, 10);
427      int i, j;
428      for (i = 1; i <= 15; ++i) for (j = 1; j <= 10; ++j)
429         A(i, j) = i + j / 1000.0;
430      DiagonalMatrix D(10);
431      D << 0.2 << 0.5 << 0.1 << 0.7 << 0.8 << 0.3 << 0.4 << 0.7 << 0.9 << 0.6;
432      Matrix U = A; Matrix V = 10 - 2 * A;
433      Matrix Prod = U * D * V.t();
434
435      DiagonalMatrix D2 = D; SortDescending(D2);
436      DiagonalMatrix D1 = D; SortSV(D1, U, V); Matrix X = D1 - D2; Print(X);
437      X = Prod - U * D1 * V.t(); Clean(X,0.000000001); Print(X);
438      U = A; V = 10 - 2 * A;
439      D1 = D; SortSV(D1, U); X = D1 - D2; Print(X);
440      D1 = D; SortSV(D1, V); X = D1 - D2; Print(X);
441      X = Prod - U * D1 * V.t(); Clean(X,0.000000001); Print(X);
442
443      D2 = D; SortAscending(D2);
444      U = A; V = 10 - 2 * A;
445      D1 = D; SortSV(D1, U, V, true); X = D1 - D2; Print(X);
446      X = Prod - U * D1 * V.t(); Clean(X,0.000000001); Print(X);
447      U = A; V = 10 - 2 * A;
448      D1 = D; SortSV(D1, U, true); X = D1 - D2; Print(X);
449      D1 = D; SortSV(D1, V, true); X = D1 - D2; Print(X);
450      X = Prod - U * D1 * V.t(); Clean(X,0.000000001); Print(X);
451   }
452
453   {
454      Tracer et1("Stage 9");
455      // Tom William's example
456      Matrix A(10,10);
457      Matrix U;
458      Matrix V;
459      DiagonalMatrix Sigma;
460      Real myVals[] =
461      {
462         1,    1,    1,    1,    1,    1,    1,    1,    1,    1,
463         1,    1,    1,    1,    1,    1,    1,    1,    1,    1,
464         1,    1,    1,    1,    1,    1,    1,    1,    1,    1,
465         1,    1,    1,    1,    1,    1,    1,    1,    1,    1,
466         1,    1,    1,    1,    1,    1,    1,    1,    1,    1,
467         1,    1,    1,    1,    1,    1,    1,    1,    1,    0,
468         1,    1,    1,    1,    1,    1,    1,    1,    1,    0,
469         1,    1,    1,    1,    1,    1,    1,    1,    0,    0,
470         1,    1,    1,    1,    1,    1,    1,    0,    0,    0,
471         1,    1,    1,    1,    1,    0,    0,    0,    0,    0,
472      };
473
474      A << myVals;
475      SVD(A, Sigma, U, V); CheckIsSorted(Sigma);
476      A -= U * Sigma * V.t();
477      Clean(A, 0.000000001); Print(A);
478   }
479
480
481
482}
Note: See TracBrowser for help on using the repository browser.