tmth.cpp

matrix library for linux and windos

      X(1,1) = BLU.LogDeterminant().Value()-B.LogDeterminant().Value();
      Clean(X,0.000000001); Print(X);

      UpperBandMatrix U; U << B; LowerBandMatrix L; L << B;
      DiagonalMatrix D; D << B;
      Print( Matrix(L + (U - D - B)) );



      for (i=1; i<=8; i++)  A.Column(i) << B.Column(i);
      Print(Matrix(A-B));
   }
   {
      Tracer et1("Stage 2");
      BandMatrix A(7,2,2);
      int i,j;
      for (i=1; i<=7; i++) for (j=1; j<=7; j++)
      {
         int k=i-j; if (k<0) k = -k;
         if (k==0) A(i,j)=6;
         else if (k==1) A(i,j) = -4;
         else if (k==2) A(i,j) = 1;
         A(1,1) = A(7,7) = 5;
      }
      DiagonalMatrix D(7); D = 0.0; A = A - D;
      BandLUMatrix B(A); Matrix M = A;
      ColumnVector V(6);
      V(1) = LogDeterminant(B).Value();
      V(2) = LogDeterminant(A).Value();
      V(3) = LogDeterminant(M).Value();
      V(4) = Determinant(B);
      V(5) = Determinant(A);
      V(6) = Determinant(M);
      V = V / 64 - 1; Clean(V,0.000000001); Print(V);
      ColumnVector X(7);

#ifdef ATandT
      Real a[7];
      // the previous statement causes a core dump in tmti.cpp
      // on the HP9000 - seems very strange. Possibly the exception
      // mechanism is failing to track the stack correctly. If you get
      // this problem replace by the following statement.
//      Real* a = new Real [7]; if (!a) exit(1);
      a[0]=1; a[1]=2; a[2]=3; a[3]=4; a[4]=5; a[5]=6; a[6]=7;
#else
      Real a[] = {1,2,3,4,5,6,7};
#endif
      X << a;
// include these if you are using the previous dynamic definition of a
//#ifdef ATandT
//      delete [] a;
//#endif
      M = (M.i()*X).t() - (B.i()*X).t() * 2.0 + (A.i()*X).t();
      Clean(M,0.000000001); Print(M);


      BandMatrix P(80,2,5); ColumnVector CX(80);
      for (i=1; i<=80; i++) for (j=1; j<=80; j++)
      { int d = i-j; if (d<=2 && d>=-5) P(i,j) = i + j/100.0; }
      for (i=1; i<=80; i++)  CX(i) = i*100.0;
      Matrix MP = P;
      ColumnVector V1 = P.i() * CX; ColumnVector V2 = MP.i() * CX;
      V = V1 - V2; Clean(V,0.000000001); Print(V);

      V1 = P * V1; V2 = MP * V2; V = V1 - V2; Clean(V,0.000000001); Print(V);
      RowVector XX(1);
      XX = LogDeterminant(P).Value() / LogDeterminant(MP).Value() - 1.0;
      Clean(XX,0.000000001); Print(XX);

      LowerBandMatrix LP(80,5);
      for (i=1; i<=80; i++) for (j=1; j<=80; j++)
      { int d = i-j; if (d<=5 && d>=0) LP(i,j) = i + j/100.0; }
      MP = LP;
      XX.ReSize(4);
      XX(1) = LogDeterminant(LP).Value();
      XX(2) = LogDeterminant(MP).Value();
      V1 = LP.i() * CX; V2 = MP.i() * CX;
      V = V1 - V2; Clean(V,0.000000001); Print(V);

      UpperBandMatrix UP(80,4);
      for (i=1; i<=80; i++) for (j=1; j<=80; j++)
      { int d = i-j; if (d<=0 && d>=-4) UP(i,j) = i + j/100.0; }
      MP = UP;
      XX(3) = LogDeterminant(UP).Value();
      XX(4) = LogDeterminant(MP).Value();
      V1 = UP.i() * CX; V2 = MP.i() * CX;
      V = V1 - V2; Clean(V,0.000000001); Print(V);
      XX = XX / SumAbsoluteValue(XX) - .25; Clean(XX,0.000000001); Print(XX);
   }

   {
      Tracer et1("Stage 3");
      SymmetricBandMatrix SA(8,5);
      int i,j;
      for (i=1; i<=8; i++) for (j=1; j<=8; j++)
         if (i-j<=5 && 0<=i-j) SA(i,j) =i + j/128.0;
      DiagonalMatrix D(8); D = 10; SA = SA + D;

      Matrix MA1(8,8); Matrix MA2(8,8);
      for (i=1; i<=8; i++)
         { MA1.Column(i) << SA.Column(i); MA2.Row(i) << SA.Row(i); }
      Print(Matrix(MA1-MA2));

      D = 10; SA = SA.t() + D; MA1 = MA1 + D;
      Print(Matrix(MA1-SA));

      UpperBandMatrix UB(8,3); LowerBandMatrix LB(8,4);
      D << SA; UB << SA; LB << SA;
      SA = SA * 5.0; D = D * 5.0; LB = LB * 5.0; UB = UB * 5.0;
      BandMatrix B = LB - D + UB - SA; Print(Matrix(B));

      SymmetricBandMatrix A(7,2); A = 100.0;
      for (i=1; i<=7; i++) for (j=1; j<=7; j++)
      {
         int k=i-j;
         if (k==0) A(i,j)=6;
         else if (k==1) A(i,j) = -4;
         else if (k==2) A(i,j) = 1;
         A(1,1) = A(7,7) = 5;
      }
      BandLUMatrix C(A); Matrix M = A;
      ColumnVector X(8);
      X(1) = LogDeterminant(C).Value() - 64;
      X(2) = LogDeterminant(A).Value() - 64;
      X(3) = LogDeterminant(M).Value() - 64;
      X(4) = SumSquare(M) - SumSquare(A);
      X(5) = SumAbsoluteValue(M) - SumAbsoluteValue(A);
      X(6) = MaximumAbsoluteValue(M) - MaximumAbsoluteValue(A);
      X(7) = Trace(M) - Trace(A);
      X(8) = Sum(M) - Sum(A);
      Clean(X,0.000000001); Print(X);

#ifdef ATandT
      Real a[7]; a[0]=1; a[1]=2; a[2]=3; a[3]=4; a[4]=5; a[5]=6; a[6]=7;
#else
      Real a[] = {1,2,3,4,5,6,7};
#endif
      X.ReSize(7);
      X << a;
      X = M.i()*X - C.i()*X * 2 + A.i()*X;
      Clean(X,0.000000001); Print(X);


      LB << A; UB << A; D << A;
      BandMatrix XA = LB + (UB - D);
      Print(Matrix(XA - A));

      for (i=1; i<=7; i++) for (j=1; j<=7; j++)
      {
         int k=i-j;
         if (k==0) A(i,j)=6;
         else if (k==1) A(i,j) = -4;
         else if (k==2) A(i,j) = 1;
         A(1,1) = A(7,7) = 5;
      }
      D = 1;

      M = LB.i() * LB - D; Clean(M,0.000000001); Print(M);
      M = UB.i() * UB - D; Clean(M,0.000000001); Print(M);
      M = XA.i() * XA - D; Clean(M,0.000000001); Print(M);
      Matrix MUB = UB; Matrix MLB = LB;
      M = LB.i() * UB - LB.i() * MUB; Clean(M,0.000000001); Print(M);
      M = UB.i() * LB - UB.i() * MLB; Clean(M,0.000000001); Print(M);
      M = LB.i() * UB - LB.i() * Matrix(UB); Clean(M,0.000000001); Print(M);
      M = UB.i() * LB - UB.i() * Matrix(LB); Clean(M,0.000000001); Print(M);
   }

   {
      // some tests about adding and subtracting band matrices of different
      // sizes - check bandwidth of results
      Tracer et1("Stage 4");

      BandFunctions(9, 3, 9, 3);   // equal
      BandFunctions(4, 7, 4, 7);   // equal
      BandFunctions(9, 3, 5, 8);   // neither < or >
      BandFunctions(5, 8, 9, 3);   // neither < or >
      BandFunctions(9, 8, 5, 3);   // >
      BandFunctions(3, 5, 8, 9);   // <

      LowerBandFunctions(9, 9);    // equal
      LowerBandFunctions(4, 4);    // equal
      LowerBandFunctions(9, 5);    // >
      LowerBandFunctions(3, 8);    // <

      UpperBandFunctions(3, 3);    // equal
      UpperBandFunctions(7, 7);    // equal
      UpperBandFunctions(8, 3);    // >
      UpperBandFunctions(5, 9);    // <

      SymmetricBandFunctions(9, 9);   // equal
      SymmetricBandFunctions(4, 4);   // equal
      SymmetricBandFunctions(9, 5);   // >
      SymmetricBandFunctions(3, 8);   // <

      DiagonalMatrix D(6); D << 2 << 3 << 4.5 << 1.25 << 9.5 << -5;
      BandMatrix BD = D;
      UpperBandMatrix UBD; UBD = D;
      LowerBandMatrix LBD; LBD = D;
      SymmetricBandMatrix SBD = D;
      Matrix X = BD - D; Print(X); X = UBD - D; Print(X);
      X = LBD - D; Print(X); X = SBD - D; Print(X);
      Matrix Test(9,2);
      Test(1,1) =  BD.BandWidth().Lower(); Test(1,2) =  BD.BandWidth().Upper();
      Test(2,1) = UBD.BandWidth().Lower(); Test(2,2) = UBD.BandWidth().Upper();
      Test(3,1) = LBD.BandWidth().Lower(); Test(3,2) = LBD.BandWidth().Upper();
      Test(4,1) = SBD.BandWidth().Lower(); Test(4,2) = SBD.BandWidth().Upper();

      IdentityMatrix I(10); I *= 5;
      BD = I; UBD = I; LBD = I; SBD = I;
      X = BD - I; Print(X); X = UBD - I; Print(X);
      X = LBD - I; Print(X); X = SBD - I; Print(X);
      Test(5,1) =  BD.BandWidth().Lower(); Test(5,2) =  BD.BandWidth().Upper();
      Test(6,1) = UBD.BandWidth().Lower(); Test(6,2) = UBD.BandWidth().Upper();
      Test(7,1) = LBD.BandWidth().Lower(); Test(7,2) = LBD.BandWidth().Upper();
      Test(8,1) = SBD.BandWidth().Lower(); Test(8,2) = SBD.BandWidth().Upper();

      RowVector RV = D.AsRow(); I.ReSize(6); BandMatrix BI = I; I = 1;
      BD =  RV.AsDiagonal() +  BI; X =  BD - D - I; Print(X);
      Test(9,1) =  BD.BandWidth().Lower(); Test(9,2) =  BD.BandWidth().Upper();

      Print(Test);
   }

   {
      // various element functions
      Tracer et1("Stage 5");

      int i, j;
      Matrix Count(1, 1); Count = 0;  // for counting errors
      Matrix M(20,30);
      for (i = 1; i <= 20; ++i) for (j = 1; j <= 30; ++j)
         M(i, j) = 100 * i + j;
      const Matrix CM = M;
      for (i = 1; i <= 20; ++i) for (j = 1; j <= 30; ++j)
      {
         if (M(i, j) != CM(i, j)) ++Count(1,1);
         if (M(i, j) != CM.element(i-1, j-1)) ++Count(1,1);
         if (M(i, j) != M.element(i-1, j-1)) ++Count(1,1);
      }

      UpperTriangularMatrix U(20);
      for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
         U(i, j) = 100 * i + j;
      const UpperTriangularMatrix CU = U;
      for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
      {
         if (U(i, j) != CU(i, j)) ++Count(1,1);
         if (U(i, j) != CU.element(i-1, j-1)) ++Count(1,1);
         if (U(i, j) != U.element(i-1, j-1)) ++Count(1,1);
      }

      LowerTriangularMatrix L(20);
      for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
         L(i, j) = 100 * i + j;
      const LowerTriangularMatrix CL = L;
      for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
      {
         if (L(i, j) != CL(i, j)) ++Count(1,1);
         if (L(i, j) != CL.element(i-1, j-1)) ++Count(1,1);
         if (L(i, j) != L.element(i-1, j-1)) ++Count(1,1);
      }

      SymmetricMatrix S(20);
      for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
         S(i, j) = 100 * i + j;
      const SymmetricMatrix CS = S;
      for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
      {
         if (S(i, j) != CS(i, j)) ++Count(1,1);
         if (S(i, j) != CS.element(i-1, j-1)) ++Count(1,1);
         if (S(i, j) != S.element(i-1, j-1)) ++Count(1,1);
         if (S(i, j) != S(j, i)) ++Count(1,1);
         if (S(i, j) != CS(i, j)) ++Count(1,1);
         if (S(i, j) != CS.element(i-1, j-1)) ++Count(1,1);
         if (S(i, j) != S.element(i-1, j-1)) ++Count(1,1);
      }

      DiagonalMatrix D(20);
      for (i = 1; i <= 20; ++i) D(i) = 100 * i + i * i;
      const DiagonalMatrix CD = D;
      for (i = 1; i <= 20; ++i)
      {
         if (D(i, i) != CD(i, i)) ++Count(1,1);
         if (D(i, i) != CD.element(i-1, i-1)) ++Count(1,1);
         if (D(i, i) != D.element(i-1, i-1)) ++Count(1,1);
         if (D(i, i) != D(i)) ++Count(1,1);
         if (D(i) != CD(i)) ++Count(1,1);
         if (D(i) != CD.element(i-1)) ++Count(1,1);
         if (D(i) != D.element(i-1)) ++Count(1,1);
      }

      RowVector R(20);
      for (i = 1; i <= 20; ++i) R(i) = 100 * i + i * i;
      const RowVector CR = R;
      for (i = 1; i <= 20; ++i)
      {
         if (R(i) != CR(i)) ++Count(1,1);
         if (R(i) != CR.element(i-1)) ++Count(1,1);
         if (R(i) != R.element(i-1)) ++Count(1,1);
      }

      ColumnVector C(20);
      for (i = 1; i <= 20; ++i) C(i) = 100 * i + i * i;
      const ColumnVector CC = C;
      for (i = 1; i <= 20; ++i)
      {
         if (C(i) != CC(i)) ++Count(1,1);
         if (C(i) != CC.element(i-1)) ++Count(1,1);
         if (C(i) != C.element(i-1)) ++Count(1,1);
      }

      Print(Count);

   }

   {
      // resize to another matrix size
      Tracer et1("Stage 6");

      Matrix A(20, 30); A = 3;
      Matrix B(3, 4);
      B.ReSize(A); B = 6; B -= 2 * A; Print(B);

      A.ReSize(25,25); A = 12;

      UpperTriangularMatrix U(5);
      U.ReSize(A); U = 12; U << (U - A); Print(U);

      LowerTriangularMatrix L(5);
      L.ReSize(U); L = 12; L << (L - A); Print(L);

      DiagonalMatrix D(5);
      D.ReSize(U); D = 12; D << (D - A); Print(D);

      SymmetricMatrix S(5);
      S.ReSize(U); S = 12; S << (S - A); Print(S);

      IdentityMatrix I(5);
      I.ReSize(U); I = 12; D << (I - A); Print(D);

      A.ReSize(10, 1); A = 17;
      ColumnVector C(5); C.ReSize(A); C = 17; C -= A; Print(C);

      A.ReSize(1, 10); A = 15;
      RowVector R(5); R.ReSize(A); R = 15; R -= A; Print(R);

   }

   {
      // generic matrix and identity matrix
      Tracer et1("Stage 7");
      IdentityMatrix I(5);
      I *= 4;
      GenericMatrix GM = I;
      GM /= 2;
      DiagonalMatrix D = GM;
      Matrix A = GM + 10;
      A -= 10;
      A -= D;
      Print(A);
   }

   {
      // SP and upper and lower triangular matrices
      Tracer et1("Stage 8");
      UpperTriangularMatrix UT(4);
      UT << 3 << 7 << 3 << 9
              << 5 << 2 << 6
                   << 8 << 0
                        << 4;
      LowerTriangularMatrix LT; LT.ReSize(UT);
      LT << 2
         << 7 << 9
         << 2 << 8 << 6
         << 1 << 0 << 3 << 5;

      DiagonalMatrix D = SP(UT, LT);
      DiagonalMatrix D1(4);
      D1 << 6 << 45 << 48 << 20;
      D -= D1; Print(D);
      BandMatrix BM = SP(UT, LT);
      Matrix X = BM - D1; Print(X);
      RowVector RV(2);
      RV(1) = BM.BandWidth().Lower();
      RV(2) = BM.BandWidth().Upper();
      Print(RV);
   }

//   cout << "\nEnd of Seventeenth test\n";
}