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


kni
Author(s): Martin Günther
autogenerated on Fri Jun 7 2019 22:06:45