00001
00002
00003
00006
00007
00008
00009
00010 #include "include.h"
00011
00012 #include "newmatap.h"
00013
00014
00015 #include "tmt.h"
00016
00017 #ifdef use_namespace
00018 using namespace NEWMAT;
00019 #endif
00020
00021 static int my_max(int p, int q) { return (p > q) ? p : q; }
00022
00023 static int my_min(int p, int q) { return (p < q) ? p : q; }
00024
00025
00026 void BandFunctions(int l1, int u1, int l2, int u2)
00027 {
00028 int i, j;
00029 BandMatrix BM1(20, l1, u1); BM1 = 999999.0;
00030 for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
00031 if (i - j <= l1 && i - j >= -u1) BM1(i, j) = 100 * i + j;
00032
00033 BandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);
00034
00035 BM2.ReSize(20, l2, u2); BM2 = 777777.0;
00036 for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
00037 if (i - j <= l2 && i - j >= -u2) BM2(i, j) = (100 * i + j) * 11;
00038
00039 BandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
00040 BMM = BM1 * BM2, BMN = -BM1;
00041
00042 Matrix M1(20,20); M1 = 0;
00043 for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
00044 if (i - j <= l1 && i - j >= -u1) M1(i, j) = 100 * i + j;
00045
00046 Matrix M2(20,20); M2 = 0;
00047 for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
00048 if (i - j <= l2 && i - j >= -u2) M2(i, j) = (100 * i + j) * 11;
00049
00050 Matrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MM = M1 * M2, MN = -M1;
00051 MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
00052 Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);
00053
00054 Matrix Test(7, 2);
00055 Test(1,1) = BM1.BandWidth().Lower() - l1;
00056 Test(1,2) = BM1.BandWidth().Upper() - u1;
00057 Test(2,1) = BM2.BandWidth().Lower() - l2;
00058 Test(2,2) = BM2.BandWidth().Upper() - u2;
00059 Test(3,1) = BMA.BandWidth().Lower() - my_max(l1,l2);
00060 Test(3,2) = BMA.BandWidth().Upper() - my_max(u1,u2);
00061 Test(4,1) = BMS.BandWidth().Lower() - my_max(l1,l2);
00062 Test(4,2) = BMS.BandWidth().Upper() - my_max(u1,u2);
00063 Test(5,1) = BMSP.BandWidth().Lower() - my_min(l1,l2);
00064 Test(5,2) = BMSP.BandWidth().Upper() - my_min(u1,u2);
00065 Test(6,1) = BMM.BandWidth().Lower() - (l1 + l2);
00066 Test(6,2) = BMM.BandWidth().Upper() - (u1 + u2);
00067 Test(7,1) = BMN.BandWidth().Lower() - l1;
00068 Test(7,2) = BMN.BandWidth().Upper() - u1;
00069 Print(Test);
00070
00071
00072 BandMatrix BM1E(20, l1, u1); BM1E = 999999.0;
00073 for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
00074 if (i - j <= l1 && i - j >= -u1) BM1E.element(i-1, j-1) = 100 * i + j;
00075 BandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
00076 for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
00077 if (i - j <= l2 && i - j >= -u2)
00078 BM2E.element(i-1, j-1) = (100 * i + j) * 11;
00079 M1 = BM1E - BM1; Print(M1);
00080 M2 = BM2E - BM2; Print(M2);
00081
00082
00083 BM1E = 444444.04; BM2E = 555555.0;
00084 const BandMatrix BM1C = BM1, BM2C = BM2;
00085 for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
00086 if (i - j <= l1 && i - j >= -u1)
00087 BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
00088 for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
00089 if (i - j <= l2 && i - j >= -u2)
00090 BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
00091 M1 = BM1E - BM1; Print(M1);
00092 M2 = BM2E - BM2; Print(M2);
00093
00094
00095 BM1E = 444444.04; BM2E = 555555.0;
00096 for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
00097 if (i - j <= l1 && i - j >= -u1) BM1E(i, j) = BM1C(i, j);
00098 for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
00099 if (i - j <= l2 && i - j >= -u2) BM2E(i, j) = BM2C(i, j);
00100 M1 = BM1E - BM1; Print(M1);
00101 M2 = BM2E - BM2; Print(M2);
00102 }
00103
00104 void LowerBandFunctions(int l1, int l2)
00105 {
00106 int i, j;
00107 LowerBandMatrix BM1(20, l1); BM1 = 999999.0;
00108 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00109 if (i - j <= l1) BM1(i, j) = 100 * i + j;
00110
00111 LowerBandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);
00112
00113 BM2.ReSize(20, l2); BM2 = 777777.0;
00114 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00115 if (i - j <= l2) BM2(i, j) = (100 * i + j) * 11;
00116
00117 LowerBandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
00118 BMM = BM1 * BM2, BMN = -BM1;
00119
00120 Matrix M1(20,20); M1 = 0;
00121 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00122 if (i - j <= l1) M1(i, j) = 100 * i + j;
00123
00124 Matrix M2(20,20); M2 = 0;
00125 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00126 if (i - j <= l2) M2(i, j) = (100 * i + j) * 11;
00127
00128 Matrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MM = M1 * M2, MN = -M1;
00129 MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
00130 Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);
00131
00132 Matrix Test(7, 2);
00133 Test(1,1) = BM1.BandWidth().Lower() - l1;
00134 Test(1,2) = BM1.BandWidth().Upper();
00135 Test(2,1) = BM2.BandWidth().Lower() - l2;
00136 Test(2,2) = BM2.BandWidth().Upper();
00137 Test(3,1) = BMA.BandWidth().Lower() - my_max(l1,l2);
00138 Test(3,2) = BMA.BandWidth().Upper();
00139 Test(4,1) = BMS.BandWidth().Lower() - my_max(l1,l2);
00140 Test(4,2) = BMS.BandWidth().Upper();
00141 Test(5,1) = BMSP.BandWidth().Lower() - my_min(l1,l2);
00142 Test(5,2) = BMSP.BandWidth().Upper();
00143 Test(6,1) = BMM.BandWidth().Lower() - (l1 + l2);
00144 Test(6,2) = BMM.BandWidth().Upper();
00145 Test(7,1) = BMN.BandWidth().Lower() - l1;
00146 Test(7,2) = BMN.BandWidth().Upper();
00147 Print(Test);
00148
00149
00150 LowerBandMatrix BM1E(20, l1); BM1E = 999999.0;
00151 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00152 if (i - j <= l1) BM1E.element(i-1, j-1) = 100 * i + j;
00153 LowerBandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
00154 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00155 if (i - j <= l2) BM2E.element(i-1, j-1) = (100 * i + j) * 11;
00156 M1 = BM1E - BM1; Print(M1);
00157 M2 = BM2E - BM2; Print(M2);
00158
00159
00160 BM1E = 444444.04; BM2E = 555555.0;
00161 const LowerBandMatrix BM1C = BM1, BM2C = BM2;
00162 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00163 if (i - j <= l1) BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
00164 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00165 if (i - j <= l2) BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
00166 M1 = BM1E - BM1; Print(M1);
00167 M2 = BM2E - BM2; Print(M2);
00168
00169
00170 BM1E = 444444.04; BM2E = 555555.0;
00171 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00172 if (i - j <= l1) BM1E(i, j) = BM1C(i, j);
00173 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00174 if (i - j <= l2) BM2E(i, j) = BM2C(i, j);
00175 M1 = BM1E - BM1; Print(M1);
00176 M2 = BM2E - BM2; Print(M2);
00177 }
00178
00179 void UpperBandFunctions(int u1, int u2)
00180 {
00181 int i, j;
00182 UpperBandMatrix BM1(20, u1); BM1 = 999999.0;
00183 for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
00184 if (i - j >= -u1) BM1(i, j) = 100 * i + j;
00185
00186 UpperBandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);
00187
00188 BM2.ReSize(20, u2); BM2 = 777777.0;
00189 for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
00190 if (i - j >= -u2) BM2(i, j) = (100 * i + j) * 11;
00191
00192 UpperBandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
00193 BMM = BM1 * BM2, BMN = -BM1;
00194
00195 Matrix M1(20,20); M1 = 0;
00196 for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
00197 if (i - j >= -u1) M1(i, j) = 100 * i + j;
00198
00199 Matrix M2(20,20); M2 = 0;
00200 for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
00201 if (i - j >= -u2) M2(i, j) = (100 * i + j) * 11;
00202
00203 Matrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MM = M1 * M2, MN = -M1;
00204 MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
00205 Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);
00206
00207 Matrix Test(7, 2);
00208 Test(1,1) = BM1.BandWidth().Lower();
00209 Test(1,2) = BM1.BandWidth().Upper() - u1;
00210 Test(2,1) = BM2.BandWidth().Lower();
00211 Test(2,2) = BM2.BandWidth().Upper() - u2;
00212 Test(3,1) = BMA.BandWidth().Lower();
00213 Test(3,2) = BMA.BandWidth().Upper() - my_max(u1,u2);
00214 Test(4,1) = BMS.BandWidth().Lower();
00215 Test(4,2) = BMS.BandWidth().Upper() - my_max(u1,u2);
00216 Test(5,1) = BMSP.BandWidth().Lower();
00217 Test(5,2) = BMSP.BandWidth().Upper() - my_min(u1,u2);
00218 Test(6,1) = BMM.BandWidth().Lower();
00219 Test(6,2) = BMM.BandWidth().Upper() - (u1 + u2);
00220 Test(7,1) = BMN.BandWidth().Lower();
00221 Test(7,2) = BMN.BandWidth().Upper() - u1;
00222 Print(Test);
00223
00224
00225 UpperBandMatrix BM1E(20, u1); BM1E = 999999.0;
00226 for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
00227 if (i - j >= -u1) BM1E.element(i-1, j-1) = 100 * i + j;
00228 UpperBandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
00229 for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
00230 if (i - j >= -u2) BM2E.element(i-1, j-1) = (100 * i + j) * 11;
00231 M1 = BM1E - BM1; Print(M1);
00232 M2 = BM2E - BM2; Print(M2);
00233
00234
00235 BM1E = 444444.04; BM2E = 555555.0;
00236 const UpperBandMatrix BM1C = BM1, BM2C = BM2;
00237 for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
00238 if (i - j >= -u1) BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
00239 for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
00240 if (i - j >= -u2) BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
00241 M1 = BM1E - BM1; Print(M1);
00242 M2 = BM2E - BM2; Print(M2);
00243
00244
00245 BM1E = 444444.04; BM2E = 555555.0;
00246 for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
00247 if (i - j >= -u1) BM1E(i, j) = BM1C(i, j);
00248 for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
00249 if (i - j >= -u2) BM2E(i, j) = BM2C(i, j);
00250 M1 = BM1E - BM1; Print(M1);
00251 M2 = BM2E - BM2; Print(M2);
00252 }
00253
00254 void SymmetricBandFunctions(int l1, int l2)
00255 {
00256 int i, j;
00257 SymmetricBandMatrix BM1(20, l1); BM1 = 999999.0;
00258 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00259 if (i - j <= l1) BM1(i, j) = 100 * i + j;
00260
00261 SymmetricBandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);
00262
00263 BM2.ReSize(20, l2); BM2 = 777777.0;
00264 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00265 if (i - j <= l2) BM2(i, j) = (100 * i + j) * 11;
00266
00267 SymmetricBandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
00268 BMN = -BM1;
00269 BandMatrix BMM = BM1 * BM2;
00270
00271 SymmetricMatrix M1(20); M1 = 0;
00272 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00273 if (i - j <= l1) M1(i, j) = 100 * i + j;
00274
00275 SymmetricMatrix M2(20); M2 = 0;
00276 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00277 if (i - j <= l2) M2(i, j) = (100 * i + j) * 11;
00278
00279 SymmetricMatrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MN = -M1;
00280 Matrix MM = M1 * M2;
00281 MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
00282 Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);
00283
00284 Matrix Test(7, 2);
00285 Test(1,1) = BM1.BandWidth().Lower() - l1;
00286 Test(1,2) = BM1.BandWidth().Upper() - l1;
00287 Test(2,1) = BM2.BandWidth().Lower() - l2;
00288 Test(2,2) = BM2.BandWidth().Upper() - l2;
00289 Test(3,1) = BMA.BandWidth().Lower() - my_max(l1,l2);
00290 Test(3,2) = BMA.BandWidth().Upper() - my_max(l1,l2);
00291 Test(4,1) = BMS.BandWidth().Lower() - my_max(l1,l2);
00292 Test(4,2) = BMS.BandWidth().Upper() - my_max(l1,l2);
00293 Test(5,1) = BMSP.BandWidth().Lower() - my_min(l1,l2);
00294 Test(5,2) = BMSP.BandWidth().Upper() - my_min(l1,l2);
00295 Test(6,1) = BMM.BandWidth().Lower() - (l1 + l2);
00296 Test(6,2) = BMM.BandWidth().Upper() - (l1 + l2);
00297 Test(7,1) = BMN.BandWidth().Lower() - l1;
00298 Test(7,2) = BMN.BandWidth().Upper() - l1;
00299 Print(Test);
00300
00301
00302 SymmetricBandMatrix BM1E(20, l1); BM1E = 999999.0;
00303 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00304 if (i - j <= l1) BM1E.element(i-1, j-1) = 100 * i + j;
00305 SymmetricBandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
00306 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00307 if (i - j <= l2) BM2E.element(i-1, j-1) = (100 * i + j) * 11;
00308 M1 = BM1E - BM1; Print(M1);
00309 M2 = BM2E - BM2; Print(M2);
00310
00311
00312 BM1E = 999999.0;
00313 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00314 if (i - j <= l1) BM1E.element(j-1, i-1) = 100 * i + j;
00315 BM2E = 777777.0;
00316 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00317 if (i - j <= l2) BM2E.element(j-1, i-1) = (100 * i + j) * 11;
00318 M1 = BM1E - BM1; Print(M1);
00319 M2 = BM2E - BM2; Print(M2);
00320
00321
00322 BM1E = 444444.04; BM2E = 555555.0;
00323 const SymmetricBandMatrix BM1C = BM1, BM2C = BM2;
00324 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00325 if (i - j <= l1) BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
00326 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00327 if (i - j <= l2) BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
00328 M1 = BM1E - BM1; Print(M1);
00329 M2 = BM2E - BM2; Print(M2);
00330
00331
00332 BM1E = 444444.0; BM2E = 555555.0;
00333 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00334 if (i - j <= l1) BM1E.element(j-1, i-1) = BM1C.element(j-1, i-1);
00335 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00336 if (i - j <= l2) BM2E.element(j-1, i-1) = BM2C.element(j-1, i-1);
00337 M1 = BM1E - BM1; Print(M1);
00338 M2 = BM2E - BM2; Print(M2);
00339
00340
00341 BM1E = 444444.0; BM2E = 555555.0;
00342 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00343 if (i - j <= l1) BM1E(i, j) = BM1C(i, j);
00344 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00345 if (i - j <= l2) BM2E(i, j) = BM2C(i, j);
00346 M1 = BM1E - BM1; Print(M1);
00347 M2 = BM2E - BM2; Print(M2);
00348
00349
00350 BM1E = 444444.0; BM2E = 555555.0;
00351 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00352 if (i - j <= l1) BM1E(j, i) = BM1C(j, i);
00353 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00354 if (i - j <= l2) BM2E(j, i) = BM2C(j, i);
00355 M1 = BM1E - BM1; Print(M1);
00356 M2 = BM2E - BM2; Print(M2);
00357
00358
00359 BM1E = 444444.0; BM2E = 555555.0;
00360 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00361 if (i - j <= l1) BM1E(j, i) = BM1C(i, j);
00362 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00363 if (i - j <= l2) BM2E(j, i) = BM2C(i, j);
00364 M1 = BM1E - BM1; Print(M1);
00365 M2 = BM2E - BM2; Print(M2);
00366
00367 }
00368
00369
00370
00371 void trymath()
00372 {
00373
00374
00375 Tracer et("Seventeenth test of Matrix package");
00376 Tracer::PrintTrace();
00377
00378
00379 {
00380 Tracer et1("Stage 1");
00381 int i, j;
00382 BandMatrix B(8,3,1);
00383 for (i=1; i<=8; i++) for (j=-3; j<=1; j++)
00384 { int k = i+j; if (k>0 && k<=8) B(i,k) = i + k/64.0; }
00385
00386 IdentityMatrix I(8); BandMatrix B1; B1 = I;
00387 UpperTriangularMatrix UT = I; Print(Matrix(B1-UT));
00388 Print(Matrix(B * B - B * 2 + I - (B - I) * (B - I)));
00389 Matrix A = B; BandMatrix C; C = B;
00390 Print(Matrix(B * A - C * 2 + I - (A - I) * (B - I)));
00391
00392 C.ReSize(8,2,2); C = 0.0; C.Inject(B);
00393 Matrix X = A.t();
00394 B1.ReSize(8,2,2); B1.Inject(X); Print(Matrix(C.t()-B1));
00395
00396 Matrix BI = B.i(); A = A.i()-BI; Clean(A,0.000000001); Print(A);
00397 BandLUMatrix BLU = B.t();
00398 BI = BLU.i(); A = X.i()-BI; Clean(A,0.000000001); Print(A);
00399 X.ReSize(1,1);
00400 X(1,1) = BLU.LogDeterminant().Value()-B.LogDeterminant().Value();
00401 Clean(X,0.000000001); Print(X);
00402
00403 UpperBandMatrix U; U << B; LowerBandMatrix L; L << B;
00404 DiagonalMatrix D; D << B;
00405 Print( Matrix(L + (U - D - B)) );
00406
00407
00408
00409 for (i=1; i<=8; i++) A.Column(i) << B.Column(i);
00410 Print(Matrix(A-B));
00411 }
00412 {
00413 Tracer et1("Stage 2");
00414 BandMatrix A(7,2,2);
00415 int i,j;
00416 for (i=1; i<=7; i++) for (j=1; j<=7; j++)
00417 {
00418 int k=i-j; if (k<0) k = -k;
00419 if (k==0) A(i,j)=6;
00420 else if (k==1) A(i,j) = -4;
00421 else if (k==2) A(i,j) = 1;
00422 A(1,1) = A(7,7) = 5;
00423 }
00424 DiagonalMatrix D(7); D = 0.0; A = A - D;
00425 BandLUMatrix B(A); Matrix M = A;
00426 ColumnVector V(6);
00427 V(1) = LogDeterminant(B).Value();
00428 V(2) = LogDeterminant(A).Value();
00429 V(3) = LogDeterminant(M).Value();
00430 V(4) = Determinant(B);
00431 V(5) = Determinant(A);
00432 V(6) = Determinant(M);
00433 V = V / 64 - 1; Clean(V,0.000000001); Print(V);
00434 ColumnVector X(7);
00435
00436 Real a[] = {1,2,3,4,5,6,7};
00437 X << a;
00438 M = (M.i()*X).t() - (B.i()*X).t() * 2.0 + (A.i()*X).t();
00439 Clean(M,0.000000001); Print(M);
00440
00441
00442 BandMatrix P(80,2,5); ColumnVector CX(80);
00443 for (i=1; i<=80; i++) for (j=1; j<=80; j++)
00444 { int d = i-j; if (d<=2 && d>=-5) P(i,j) = i + j/100.0; }
00445 for (i=1; i<=80; i++) CX(i) = i*100.0;
00446 Matrix MP = P;
00447 ColumnVector V1 = P.i() * CX; ColumnVector V2 = MP.i() * CX;
00448 V = V1 - V2; Clean(V,0.000000001); Print(V);
00449
00450 V1 = P * V1; V2 = MP * V2; V = V1 - V2; Clean(V,0.000000001); Print(V);
00451 RowVector XX(1);
00452 XX = LogDeterminant(P).Value() / LogDeterminant(MP).Value() - 1.0;
00453 Clean(XX,0.000000001); Print(XX);
00454
00455 LowerBandMatrix LP(80,5);
00456 for (i=1; i<=80; i++) for (j=1; j<=80; j++)
00457 { int d = i-j; if (d<=5 && d>=0) LP(i,j) = i + j/100.0; }
00458 MP = LP;
00459 XX.ReSize(4);
00460 XX(1) = LogDeterminant(LP).Value();
00461 XX(2) = LogDeterminant(MP).Value();
00462 V1 = LP.i() * CX; V2 = MP.i() * CX;
00463 V = V1 - V2; Clean(V,0.000000001); Print(V);
00464
00465 UpperBandMatrix UP(80,4);
00466 for (i=1; i<=80; i++) for (j=1; j<=80; j++)
00467 { int d = i-j; if (d<=0 && d>=-4) UP(i,j) = i + j/100.0; }
00468 MP = UP;
00469 XX(3) = LogDeterminant(UP).Value();
00470 XX(4) = LogDeterminant(MP).Value();
00471 V1 = UP.i() * CX; V2 = MP.i() * CX;
00472 V = V1 - V2; Clean(V,0.000000001); Print(V);
00473 XX = XX / SumAbsoluteValue(XX) - .25; Clean(XX,0.000000001); Print(XX);
00474 }
00475
00476 {
00477 Tracer et1("Stage 3");
00478 SymmetricBandMatrix SA(8,5);
00479 int i,j;
00480 for (i=1; i<=8; i++) for (j=1; j<=8; j++)
00481 if (i-j<=5 && 0<=i-j) SA(i,j) =i + j/128.0;
00482 DiagonalMatrix D(8); D = 10; SA = SA + D;
00483
00484 Matrix MA1(8,8); Matrix MA2(8,8);
00485 for (i=1; i<=8; i++)
00486 { MA1.Column(i) << SA.Column(i); MA2.Row(i) << SA.Row(i); }
00487 Print(Matrix(MA1-MA2));
00488
00489 D = 10; SA = SA.t() + D; MA1 = MA1 + D;
00490 Print(Matrix(MA1-SA));
00491
00492 UpperBandMatrix UB(8,3); LowerBandMatrix LB(8,4);
00493 D << SA; UB << SA; LB << SA;
00494 SA = SA * 5.0; D = D * 5.0; LB = LB * 5.0; UB = UB * 5.0;
00495 BandMatrix B = LB - D + UB - SA; Print(Matrix(B));
00496
00497 SymmetricBandMatrix A(7,2); A = 100.0;
00498 for (i=1; i<=7; i++) for (j=1; j<=7; j++)
00499 {
00500 int k=i-j;
00501 if (k==0) A(i,j)=6;
00502 else if (k==1) A(i,j) = -4;
00503 else if (k==2) A(i,j) = 1;
00504 A(1,1) = A(7,7) = 5;
00505 }
00506 BandLUMatrix C(A); Matrix M = A;
00507 ColumnVector X(8);
00508 X(1) = LogDeterminant(C).Value() - 64;
00509 X(2) = LogDeterminant(A).Value() - 64;
00510 X(3) = LogDeterminant(M).Value() - 64;
00511 X(4) = SumSquare(M) - SumSquare(A);
00512 X(5) = SumAbsoluteValue(M) - SumAbsoluteValue(A);
00513 X(6) = MaximumAbsoluteValue(M) - MaximumAbsoluteValue(A);
00514 X(7) = Trace(M) - Trace(A);
00515 X(8) = Sum(M) - Sum(A);
00516 Clean(X,0.000000001); Print(X);
00517
00518 Real a[] = {1,2,3,4,5,6,7};
00519 X.ReSize(7);
00520 X << a;
00521 X = M.i()*X - C.i()*X * 2 + A.i()*X;
00522 Clean(X,0.000000001); Print(X);
00523
00524
00525 LB << A; UB << A; D << A;
00526 BandMatrix XA = LB + (UB - D);
00527 Print(Matrix(XA - A));
00528
00529 for (i=1; i<=7; i++) for (j=1; j<=7; j++)
00530 {
00531 int k=i-j;
00532 if (k==0) A(i,j)=6;
00533 else if (k==1) A(i,j) = -4;
00534 else if (k==2) A(i,j) = 1;
00535 A(1,1) = A(7,7) = 5;
00536 }
00537 D = 1;
00538
00539 M = LB.i() * LB - D; Clean(M,0.000000001); Print(M);
00540 M = UB.i() * UB - D; Clean(M,0.000000001); Print(M);
00541 M = XA.i() * XA - D; Clean(M,0.000000001); Print(M);
00542 Matrix MUB = UB; Matrix MLB = LB;
00543 M = LB.i() * UB - LB.i() * MUB; Clean(M,0.000000001); Print(M);
00544 M = UB.i() * LB - UB.i() * MLB; Clean(M,0.000000001); Print(M);
00545 M = LB.i() * UB - LB.i() * Matrix(UB); Clean(M,0.000000001); Print(M);
00546 M = UB.i() * LB - UB.i() * Matrix(LB); Clean(M,0.000000001); Print(M);
00547 }
00548
00549 {
00550
00551
00552 Tracer et1("Stage 4");
00553
00554 BandFunctions(9, 3, 9, 3);
00555 BandFunctions(4, 7, 4, 7);
00556 BandFunctions(9, 3, 5, 8);
00557 BandFunctions(5, 8, 9, 3);
00558 BandFunctions(9, 8, 5, 3);
00559 BandFunctions(3, 5, 8, 9);
00560
00561 LowerBandFunctions(9, 9);
00562 LowerBandFunctions(4, 4);
00563 LowerBandFunctions(9, 5);
00564 LowerBandFunctions(3, 8);
00565
00566 UpperBandFunctions(3, 3);
00567 UpperBandFunctions(7, 7);
00568 UpperBandFunctions(8, 3);
00569 UpperBandFunctions(5, 9);
00570
00571 SymmetricBandFunctions(9, 9);
00572 SymmetricBandFunctions(4, 4);
00573 SymmetricBandFunctions(9, 5);
00574 SymmetricBandFunctions(3, 8);
00575
00576 DiagonalMatrix D(6); D << 2 << 3 << 4.5 << 1.25 << 9.5 << -5;
00577 BandMatrix BD = D;
00578 UpperBandMatrix UBD; UBD = D;
00579 LowerBandMatrix LBD; LBD = D;
00580 SymmetricBandMatrix SBD = D;
00581 Matrix X = BD - D; Print(X); X = UBD - D; Print(X);
00582 X = LBD - D; Print(X); X = SBD - D; Print(X);
00583 Matrix Test(9,2);
00584 Test(1,1) = BD.BandWidth().Lower(); Test(1,2) = BD.BandWidth().Upper();
00585 Test(2,1) = UBD.BandWidth().Lower(); Test(2,2) = UBD.BandWidth().Upper();
00586 Test(3,1) = LBD.BandWidth().Lower(); Test(3,2) = LBD.BandWidth().Upper();
00587 Test(4,1) = SBD.BandWidth().Lower(); Test(4,2) = SBD.BandWidth().Upper();
00588
00589 IdentityMatrix I(10); I *= 5;
00590 BD = I; UBD = I; LBD = I; SBD = I;
00591 X = BD - I; Print(X); X = UBD - I; Print(X);
00592 X = LBD - I; Print(X); X = SBD - I; Print(X);
00593 Test(5,1) = BD.BandWidth().Lower(); Test(5,2) = BD.BandWidth().Upper();
00594 Test(6,1) = UBD.BandWidth().Lower(); Test(6,2) = UBD.BandWidth().Upper();
00595 Test(7,1) = LBD.BandWidth().Lower(); Test(7,2) = LBD.BandWidth().Upper();
00596 Test(8,1) = SBD.BandWidth().Lower(); Test(8,2) = SBD.BandWidth().Upper();
00597
00598 RowVector RV = D.AsRow(); I.ReSize(6); BandMatrix BI = I; I = 1;
00599 BD = RV.AsDiagonal() + BI; X = BD - D - I; Print(X);
00600 Test(9,1) = BD.BandWidth().Lower(); Test(9,2) = BD.BandWidth().Upper();
00601
00602 Print(Test);
00603 }
00604
00605 {
00606
00607 Tracer et1("Stage 5");
00608
00609 int i, j;
00610 Matrix Count(1, 1); Count = 0;
00611 Matrix M(20,30);
00612 for (i = 1; i <= 20; ++i) for (j = 1; j <= 30; ++j)
00613 M(i, j) = 100 * i + j;
00614 const Matrix CM = M;
00615 for (i = 1; i <= 20; ++i) for (j = 1; j <= 30; ++j)
00616 {
00617 if (M(i, j) != CM(i, j)) ++Count(1,1);
00618 if (M(i, j) != CM.element(i-1, j-1)) ++Count(1,1);
00619 if (M(i, j) != M.element(i-1, j-1)) ++Count(1,1);
00620 }
00621
00622 UpperTriangularMatrix U(20);
00623 for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
00624 U(i, j) = 100 * i + j;
00625 const UpperTriangularMatrix CU = U;
00626 for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
00627 {
00628 if (U(i, j) != CU(i, j)) ++Count(1,1);
00629 if (U(i, j) != CU.element(i-1, j-1)) ++Count(1,1);
00630 if (U(i, j) != U.element(i-1, j-1)) ++Count(1,1);
00631 }
00632
00633 LowerTriangularMatrix L(20);
00634 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00635 L(i, j) = 100 * i + j;
00636 const LowerTriangularMatrix CL = L;
00637 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00638 {
00639 if (L(i, j) != CL(i, j)) ++Count(1,1);
00640 if (L(i, j) != CL.element(i-1, j-1)) ++Count(1,1);
00641 if (L(i, j) != L.element(i-1, j-1)) ++Count(1,1);
00642 }
00643
00644 SymmetricMatrix S(20);
00645 for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
00646 S(i, j) = 100 * i + j;
00647 const SymmetricMatrix CS = S;
00648 for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
00649 {
00650 if (S(i, j) != CS(i, j)) ++Count(1,1);
00651 if (S(i, j) != CS.element(i-1, j-1)) ++Count(1,1);
00652 if (S(i, j) != S.element(i-1, j-1)) ++Count(1,1);
00653 if (S(i, j) != S(j, i)) ++Count(1,1);
00654 if (S(i, j) != CS(i, j)) ++Count(1,1);
00655 if (S(i, j) != CS.element(i-1, j-1)) ++Count(1,1);
00656 if (S(i, j) != S.element(i-1, j-1)) ++Count(1,1);
00657 }
00658
00659 DiagonalMatrix D(20);
00660 for (i = 1; i <= 20; ++i) D(i) = 100 * i + i * i;
00661 const DiagonalMatrix CD = D;
00662 for (i = 1; i <= 20; ++i)
00663 {
00664 if (D(i, i) != CD(i, i)) ++Count(1,1);
00665 if (D(i, i) != CD.element(i-1, i-1)) ++Count(1,1);
00666 if (D(i, i) != D.element(i-1, i-1)) ++Count(1,1);
00667 if (D(i, i) != D(i)) ++Count(1,1);
00668 if (D(i) != CD(i)) ++Count(1,1);
00669 if (D(i) != CD.element(i-1)) ++Count(1,1);
00670 if (D(i) != D.element(i-1)) ++Count(1,1);
00671 }
00672
00673 RowVector R(20);
00674 for (i = 1; i <= 20; ++i) R(i) = 100 * i + i * i;
00675 const RowVector CR = R;
00676 for (i = 1; i <= 20; ++i)
00677 {
00678 if (R(i) != CR(i)) ++Count(1,1);
00679 if (R(i) != CR.element(i-1)) ++Count(1,1);
00680 if (R(i) != R.element(i-1)) ++Count(1,1);
00681 }
00682
00683 ColumnVector C(20);
00684 for (i = 1; i <= 20; ++i) C(i) = 100 * i + i * i;
00685 const ColumnVector CC = C;
00686 for (i = 1; i <= 20; ++i)
00687 {
00688 if (C(i) != CC(i)) ++Count(1,1);
00689 if (C(i) != CC.element(i-1)) ++Count(1,1);
00690 if (C(i) != C.element(i-1)) ++Count(1,1);
00691 }
00692
00693 Print(Count);
00694
00695 }
00696
00697 {
00698
00699 Tracer et1("Stage 6");
00700
00701 Matrix A(20, 30); A = 3;
00702 Matrix B(3, 4);
00703 B.ReSize(A); B = 6; B -= 2 * A; Print(B);
00704
00705 A.ReSize(25,25); A = 12;
00706
00707 UpperTriangularMatrix U(5);
00708 U.ReSize(A); U = 12; U << (U - A); Print(U);
00709
00710 LowerTriangularMatrix L(5);
00711 L.ReSize(U); L = 12; L << (L - A); Print(L);
00712
00713 DiagonalMatrix D(5);
00714 D.ReSize(U); D = 12; D << (D - A); Print(D);
00715
00716 SymmetricMatrix S(5);
00717 S.ReSize(U); S = 12; S << (S - A); Print(S);
00718
00719 IdentityMatrix I(5);
00720 I.ReSize(U); I = 12; D << (I - A); Print(D);
00721
00722 A.ReSize(10, 1); A = 17;
00723 ColumnVector C(5); C.ReSize(A); C = 17; C -= A; Print(C);
00724
00725 A.ReSize(1, 10); A = 15;
00726 RowVector R(5); R.ReSize(A); R = 15; R -= A; Print(R);
00727
00728 }
00729
00730 {
00731
00732 Tracer et1("Stage 7");
00733 IdentityMatrix I(5);
00734 I *= 4;
00735 GenericMatrix GM = I;
00736 GM /= 2;
00737 DiagonalMatrix D = GM;
00738 Matrix A = GM + 10;
00739 A -= 10;
00740 A -= D;
00741 Print(A);
00742 }
00743
00744 {
00745
00746 Tracer et1("Stage 8");
00747 UpperTriangularMatrix UT(4);
00748 UT << 3 << 7 << 3 << 9
00749 << 5 << 2 << 6
00750 << 8 << 0
00751 << 4;
00752 LowerTriangularMatrix LT; LT.ReSize(UT);
00753 LT << 2
00754 << 7 << 9
00755 << 2 << 8 << 6
00756 << 1 << 0 << 3 << 5;
00757
00758 DiagonalMatrix D = SP(UT, LT);
00759 DiagonalMatrix D1(4);
00760 D1 << 6 << 45 << 48 << 20;
00761 D -= D1; Print(D);
00762 BandMatrix BM = SP(UT, LT);
00763 Matrix X = BM - D1; Print(X);
00764 RowVector RV(2);
00765 RV(1) = BM.BandWidth().Lower();
00766 RV(2) = BM.BandWidth().Upper();
00767 Print(RV);
00768 }
00769
00770 {
00771
00772 Tracer et1("Stage 9");
00773 MultWithCarry MCW;
00774 int i, j;
00775
00776 IdentityMatrix I(8);
00777 Matrix X = I;
00778 Matrix Y = Helmert_transpose(X);
00779 Matrix H = Helmert(9); H -= Y.t(); Clean(H,0.000000001); Print(H);
00780 Matrix Z = Helmert(Y) - I;
00781 Clean(Z,0.000000001); Print(Z);
00782
00783 Matrix A(9, 8);
00784 for (i = 1; i <= 9; ++i) for (j = 1; j <= 8; ++j)
00785 A(i, j) = Helmert_transpose(X.column(j), i);
00786 A -= Y; Clean(A,0.000000001); Print(A);
00787
00788 X = I;
00789 Y = Helmert_transpose(X, true);
00790 H = Helmert(8, true); H -= Y.t(); Clean(H,0.000000001); Print(H);
00791 Z = Helmert(Y, true) - I;
00792 Clean(Z,0.000000001); Print(Z);
00793
00794 A.resize(8, 8);
00795 for (i = 1; i <= 8; ++i) for (j = 1; j <= 8; ++j)
00796 A(i, j) = Helmert_transpose(X.column(j), i, true);
00797 A -= Y; Clean(A,0.000000001); Print(A);
00798
00799
00800
00801 I.ReSize(9);
00802 X = I;
00803 Y = Helmert(X, true);
00804 H = Helmert(9, true); H -= Y; Clean(H,0.000000001); Print(H);
00805 Z = Helmert_transpose(Y, true) - I;
00806 Clean(Z,0.000000001); Print(Z);
00807
00808 A.ReSize(9, 9);
00809 for (i = 1; i <= 9; ++i) A.Column(i) = Helmert(9, i, true);
00810 A -= Y; Clean(A,0.000000001); Print(A);
00811
00812 Y = Helmert(X);
00813 A.ReSize(8, 9);
00814 for (i = 1; i <= 9; ++i) A.Column(i) = Helmert(9, i);
00815 A -= Y; Clean(A,0.000000001); Print(A);
00816
00817 ColumnVector Twos(100); Twos = 2;
00818 ColumnVector CV = Helmert(Twos); Clean(CV,0.000000001); Print(CV);
00819
00820 X.resize(25,30);
00821 FillWithValues(MCW, X);
00822 Y = Helmert(X);
00823 Z = Helmert(X,true).rows(1,24) - Y;
00824 Clean(Z,0.000000001); Print(Z);
00825 Z = Helmert(X,true).row(25) - X.sum_columns() / 5.0;
00826 Clean(Z,0.000000001); Print(Z);
00827
00828 I.resize(15);
00829 X = I;
00830 Z = Helmert_transpose(X, true) - Helmert(X, true).t();
00831 Clean(Z,0.000000001); Print(Z);
00832 I.resize(14); Y = I;
00833 Z = Helmert(X) - Helmert_transpose(Y).t();
00834 Clean(Z,0.000000001); Print(Z);
00835
00836
00837
00838 }
00839
00840
00841
00842
00843
00844
00845
00846 }
00847
00848