[4565] | 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 |
---|
| 13 | using namespace NEWMAT; |
---|
| 14 | #endif |
---|
| 15 | |
---|
| 16 | // check D is sorted |
---|
| 17 | void 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 | |
---|
| 26 | void 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 | } |
---|