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

rl_agent
Author(s): Todd Hester
autogenerated on Thu Jun 6 2019 22:00:14