tmth.cpp

matrix library for linux and windos

//#define WANT_STREAM

#include "include.h"

#include "newmatap.h"
//#include "newmatio.h"

#include "tmt.h"

#ifdef use_namespace
using namespace NEWMAT;
#endif

static int my_max(int p, int q) { return (p > q) ? p : q; }

static int my_min(int p, int q) { return (p < q) ? p : q; }


void BandFunctions(int l1, int u1, int l2, int u2)
{
   int i, j;
   BandMatrix BM1(20, l1, u1); BM1 = 999999.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
      if (i - j <= l1 && i - j >= -u1) BM1(i, j) = 100 * i + j;

   BandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);

   BM2.ReSize(20, l2, u2); BM2 = 777777.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
      if (i - j <= l2 && i - j >= -u2) BM2(i, j) = (100 * i + j) * 11;

   BandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
      BMM = BM1 * BM2, BMN = -BM1;

   Matrix M1(20,20); M1 = 0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
      if (i - j <= l1 && i - j >= -u1) M1(i, j) = 100 * i + j;

   Matrix M2(20,20); M2 = 0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
      if (i - j <= l2 && i - j >= -u2) M2(i, j) = (100 * i + j) * 11;

   Matrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MM = M1 * M2, MN = -M1;
   MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
   Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);

   Matrix Test(7, 2);
   Test(1,1) = BM1.BandWidth().Lower() - l1;
   Test(1,2) = BM1.BandWidth().Upper() - u1;
   Test(2,1) = BM2.BandWidth().Lower() - l2;
   Test(2,2) = BM2.BandWidth().Upper() - u2;
   Test(3,1) = BMA.BandWidth().Lower() - my_max(l1,l2);
   Test(3,2) = BMA.BandWidth().Upper() - my_max(u1,u2);
   Test(4,1) = BMS.BandWidth().Lower() - my_max(l1,l2);
   Test(4,2) = BMS.BandWidth().Upper() - my_max(u1,u2);
   Test(5,1) = BMSP.BandWidth().Lower() - my_min(l1,l2);
   Test(5,2) = BMSP.BandWidth().Upper() - my_min(u1,u2);
   Test(6,1) = BMM.BandWidth().Lower() - (l1 + l2);
   Test(6,2) = BMM.BandWidth().Upper() - (u1 + u2);
   Test(7,1) = BMN.BandWidth().Lower() - l1;
   Test(7,2) = BMN.BandWidth().Upper() - u1;
   Print(Test);

   // test element function
   BandMatrix BM1E(20, l1, u1); BM1E = 999999.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
      if (i - j <= l1 && i - j >= -u1) BM1E.element(i-1, j-1) = 100 * i + j;
   BandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
      if (i - j <= l2 && i - j >= -u2)
         BM2E.element(i-1, j-1) = (100 * i + j) * 11;
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);

   // test element function with constant
   BM1E = 444444.04; BM2E = 555555.0;
   const BandMatrix BM1C = BM1, BM2C = BM2;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
      if (i - j <= l1 && i - j >= -u1)
         BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
   for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
      if (i - j <= l2 && i - j >= -u2)
         BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);

   // test subscript function with constant
   BM1E = 444444.04; BM2E = 555555.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
      if (i - j <= l1 && i - j >= -u1) BM1E(i, j) = BM1C(i, j);
   for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
      if (i - j <= l2 && i - j >= -u2) BM2E(i, j) = BM2C(i, j);
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);
}

void LowerBandFunctions(int l1, int l2)
{
   int i, j;
   LowerBandMatrix BM1(20, l1); BM1 = 999999.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l1) BM1(i, j) = 100 * i + j;

   LowerBandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);

   BM2.ReSize(20, l2); BM2 = 777777.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l2) BM2(i, j) = (100 * i + j) * 11;

   LowerBandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
      BMM = BM1 * BM2, BMN = -BM1;

   Matrix M1(20,20); M1 = 0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l1) M1(i, j) = 100 * i + j;

   Matrix M2(20,20); M2 = 0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l2) M2(i, j) = (100 * i + j) * 11;

   Matrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MM = M1 * M2, MN = -M1;
   MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
   Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);

   Matrix Test(7, 2);
   Test(1,1) = BM1.BandWidth().Lower() - l1;
   Test(1,2) = BM1.BandWidth().Upper();
   Test(2,1) = BM2.BandWidth().Lower() - l2;
   Test(2,2) = BM2.BandWidth().Upper();
   Test(3,1) = BMA.BandWidth().Lower() - my_max(l1,l2);
   Test(3,2) = BMA.BandWidth().Upper();
   Test(4,1) = BMS.BandWidth().Lower() - my_max(l1,l2);
   Test(4,2) = BMS.BandWidth().Upper();
   Test(5,1) = BMSP.BandWidth().Lower() - my_min(l1,l2);
   Test(5,2) = BMSP.BandWidth().Upper();
   Test(6,1) = BMM.BandWidth().Lower() - (l1 + l2);
   Test(6,2) = BMM.BandWidth().Upper();
   Test(7,1) = BMN.BandWidth().Lower() - l1;
   Test(7,2) = BMN.BandWidth().Upper();
   Print(Test);

   // test element function
   LowerBandMatrix BM1E(20, l1); BM1E = 999999.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l1) BM1E.element(i-1, j-1) = 100 * i + j;
   LowerBandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l2) BM2E.element(i-1, j-1) = (100 * i + j) * 11;
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);

   // test element function with constant
   BM1E = 444444.04; BM2E = 555555.0;
   const LowerBandMatrix BM1C = BM1, BM2C = BM2;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l1) BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l2) BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);

   // test subscript function with constant
   BM1E = 444444.04; BM2E = 555555.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l1) BM1E(i, j) = BM1C(i, j);
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l2) BM2E(i, j) = BM2C(i, j);
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);
}

void UpperBandFunctions(int u1, int u2)
{
   int i, j;
   UpperBandMatrix BM1(20, u1); BM1 = 999999.0;
   for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
      if (i - j >= -u1) BM1(i, j) = 100 * i + j;

   UpperBandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);

   BM2.ReSize(20, u2); BM2 = 777777.0;
   for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
      if (i - j >= -u2) BM2(i, j) = (100 * i + j) * 11;

   UpperBandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
      BMM = BM1 * BM2, BMN = -BM1;

   Matrix M1(20,20); M1 = 0;
   for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
      if (i - j >= -u1) M1(i, j) = 100 * i + j;

   Matrix M2(20,20); M2 = 0;
   for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
      if (i - j >= -u2) M2(i, j) = (100 * i + j) * 11;

   Matrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MM = M1 * M2, MN = -M1;
   MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
   Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);

   Matrix Test(7, 2);
   Test(1,1) = BM1.BandWidth().Lower();
   Test(1,2) = BM1.BandWidth().Upper() - u1;
   Test(2,1) = BM2.BandWidth().Lower();
   Test(2,2) = BM2.BandWidth().Upper() - u2;
   Test(3,1) = BMA.BandWidth().Lower();
   Test(3,2) = BMA.BandWidth().Upper() - my_max(u1,u2);
   Test(4,1) = BMS.BandWidth().Lower();
   Test(4,2) = BMS.BandWidth().Upper() - my_max(u1,u2);
   Test(5,1) = BMSP.BandWidth().Lower();
   Test(5,2) = BMSP.BandWidth().Upper() - my_min(u1,u2);
   Test(6,1) = BMM.BandWidth().Lower();
   Test(6,2) = BMM.BandWidth().Upper() - (u1 + u2);
   Test(7,1) = BMN.BandWidth().Lower();
   Test(7,2) = BMN.BandWidth().Upper() - u1;
   Print(Test);

   // test element function
   UpperBandMatrix BM1E(20, u1); BM1E = 999999.0;
   for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
      if (i - j >= -u1) BM1E.element(i-1, j-1) = 100 * i + j;
   UpperBandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
   for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
      if (i - j >= -u2) BM2E.element(i-1, j-1) = (100 * i + j) * 11;
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);

   // test element function with constant
   BM1E = 444444.04; BM2E = 555555.0;
   const UpperBandMatrix BM1C = BM1, BM2C = BM2;
   for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
      if (i - j >= -u1) BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
   for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
      if (i - j >= -u2) BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);

   // test subscript function with constant
   BM1E = 444444.04; BM2E = 555555.0;
   for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
      if (i - j >= -u1) BM1E(i, j) = BM1C(i, j);
   for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
      if (i - j >= -u2) BM2E(i, j) = BM2C(i, j);
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);
}

void SymmetricBandFunctions(int l1, int l2)
{
   int i, j;
   SymmetricBandMatrix BM1(20, l1); BM1 = 999999.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l1) BM1(i, j) = 100 * i + j;

   SymmetricBandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);

   BM2.ReSize(20, l2); BM2 = 777777.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l2) BM2(i, j) = (100 * i + j) * 11;

   SymmetricBandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
      BMN = -BM1;
   BandMatrix BMM = BM1 * BM2;

   SymmetricMatrix M1(20); M1 = 0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l1) M1(i, j) = 100 * i + j;

   SymmetricMatrix M2(20); M2 = 0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l2) M2(i, j) = (100 * i + j) * 11;

   SymmetricMatrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MN = -M1;
   Matrix MM = M1 * M2;
   MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
   Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);

   Matrix Test(7, 2);
   Test(1,1) = BM1.BandWidth().Lower() - l1;
   Test(1,2) = BM1.BandWidth().Upper() - l1;
   Test(2,1) = BM2.BandWidth().Lower() - l2;
   Test(2,2) = BM2.BandWidth().Upper() - l2;
   Test(3,1) = BMA.BandWidth().Lower() - my_max(l1,l2);
   Test(3,2) = BMA.BandWidth().Upper() - my_max(l1,l2);
   Test(4,1) = BMS.BandWidth().Lower() - my_max(l1,l2);
   Test(4,2) = BMS.BandWidth().Upper() - my_max(l1,l2);
   Test(5,1) = BMSP.BandWidth().Lower() - my_min(l1,l2);
   Test(5,2) = BMSP.BandWidth().Upper() - my_min(l1,l2);
   Test(6,1) = BMM.BandWidth().Lower() - (l1 + l2);
   Test(6,2) = BMM.BandWidth().Upper() - (l1 + l2);
   Test(7,1) = BMN.BandWidth().Lower() - l1;
   Test(7,2) = BMN.BandWidth().Upper() - l1;
   Print(Test);

   // test element function
   SymmetricBandMatrix BM1E(20, l1); BM1E = 999999.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l1) BM1E.element(i-1, j-1) = 100 * i + j;
   SymmetricBandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l2) BM2E.element(i-1, j-1) = (100 * i + j) * 11;
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);

   // reverse subscripts
   BM1E = 999999.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l1) BM1E.element(j-1, i-1) = 100 * i + j;
   BM2E = 777777.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l2) BM2E.element(j-1, i-1) = (100 * i + j) * 11;
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);

   // test element function with constant
   BM1E = 444444.04; BM2E = 555555.0;
   const SymmetricBandMatrix BM1C = BM1, BM2C = BM2;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l1) BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l2) BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);

   // reverse subscripts
   BM1E = 444444.0; BM2E = 555555.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l1) BM1E.element(j-1, i-1) = BM1C.element(j-1, i-1);
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l2) BM2E.element(j-1, i-1) = BM2C.element(j-1, i-1);
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);

   // test subscript function with constant
   BM1E = 444444.0; BM2E = 555555.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l1) BM1E(i, j) = BM1C(i, j);
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l2) BM2E(i, j) = BM2C(i, j);
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);

   // reverse subscripts
   BM1E = 444444.0; BM2E = 555555.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l1) BM1E(j, i) = BM1C(j, i);
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l2) BM2E(j, i) = BM2C(j, i);
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);

   // partly reverse subscripts
   BM1E = 444444.0; BM2E = 555555.0;
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l1) BM1E(j, i) = BM1C(i, j);
   for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      if (i - j <= l2) BM2E(j, i) = BM2C(i, j);
   M1 = BM1E - BM1; Print(M1);
   M2 = BM2E - BM2; Print(M2);

}



void trymath()
{
//   cout << "\nSeventeenth test of Matrix package\n";
//   cout << "\n";
   Tracer et("Seventeenth test of Matrix package");
   Tracer::PrintTrace();


   {
      Tracer et1("Stage 1");
      int i, j;
      BandMatrix B(8,3,1);
      for (i=1; i<=8; i++) for (j=-3; j<=1; j++)
         { int k = i+j; if (k>0 && k<=8) B(i,k) = i + k/64.0; }

      IdentityMatrix I(8); BandMatrix B1; B1 = I;
      UpperTriangularMatrix UT = I; Print(Matrix(B1-UT));
      Print(Matrix(B * B - B * 2 + I - (B - I) * (B - I)));
      Matrix A = B; BandMatrix C; C = B;
      Print(Matrix(B * A - C * 2 + I - (A - I) * (B - I)));

      C.ReSize(8,2,2); C = 0.0; C.Inject(B);
      Matrix X = A.t();
      B1.ReSize(8,2,2); B1.Inject(X); Print(Matrix(C.t()-B1));

      Matrix BI = B.i(); A = A.i()-BI; Clean(A,0.000000001); Print(A);
      BandLUMatrix BLU = B.t();
      BI = BLU.i(); A = X.i()-BI; Clean(A,0.000000001); Print(A);
      X.ReSize(1,1);