tmth.cpp
Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 
00008 //#define WANT_STREAM
00009 
00010 #include "include.h"
00011 
00012 #include "newmatap.h"
00013 //#include "newmatio.h"
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    // test element function
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    // test element function with constant
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    // test subscript function with constant
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    // test element function
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    // test element function with constant
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    // test subscript function with constant
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    // test element function
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    // test element function with constant
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    // test subscript function with constant
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    // test element function
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    // reverse subscripts
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    // test element function with constant
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    // reverse subscripts
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    // test subscript function with constant
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    // reverse subscripts
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    // partly reverse subscripts
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 //   cout << "\nSeventeenth test of Matrix package\n";
00374 //   cout << "\n";
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       // some tests about adding and subtracting band matrices of different
00551       // sizes - check bandwidth of results
00552       Tracer et1("Stage 4");
00553 
00554       BandFunctions(9, 3, 9, 3);   // equal
00555       BandFunctions(4, 7, 4, 7);   // equal
00556       BandFunctions(9, 3, 5, 8);   // neither < or >
00557       BandFunctions(5, 8, 9, 3);   // neither < or >
00558       BandFunctions(9, 8, 5, 3);   // >
00559       BandFunctions(3, 5, 8, 9);   // <
00560 
00561       LowerBandFunctions(9, 9);    // equal
00562       LowerBandFunctions(4, 4);    // equal
00563       LowerBandFunctions(9, 5);    // >
00564       LowerBandFunctions(3, 8);    // <
00565 
00566       UpperBandFunctions(3, 3);    // equal
00567       UpperBandFunctions(7, 7);    // equal
00568       UpperBandFunctions(8, 3);    // >
00569       UpperBandFunctions(5, 9);    // <
00570 
00571       SymmetricBandFunctions(9, 9);   // equal
00572       SymmetricBandFunctions(4, 4);   // equal
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       // various element functions
00607       Tracer et1("Stage 5");
00608 
00609       int i, j;
00610       Matrix Count(1, 1); Count = 0;  // for counting errors
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       // resize to another matrix size
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       // generic matrix and identity matrix
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       // SP and upper and lower triangular matrices
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       // Helmert multiplies
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 //   cout << "\nEnd of Seventeenth test\n";
00846 }
00847 
00848

kni
Author(s): Martin Günther
autogenerated on Mon Aug 14 2017 02:44:13