home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / misc / volume34 / newmat06 / part07 / tmth.cxx < prev    next >
C/C++ Source or Header  |  1992-12-06  |  6KB  |  201 lines
  1.  
  2. //#define WANT_STREAM
  3.  
  4. #include "include.h"
  5.  
  6. #include "newmatap.h"
  7. //#include "newmatio.h"
  8.  
  9. void Print(const Matrix& X);
  10. void Print(const UpperTriangularMatrix& X);
  11. void Print(const DiagonalMatrix& X);
  12. void Print(const SymmetricMatrix& X);
  13. void Print(const LowerTriangularMatrix& X);
  14.  
  15. void Clean(Matrix&, Real);
  16.  
  17.  
  18.  
  19.  
  20. void trymath()
  21. {
  22. //   cout << "\nSeventeenth test of Matrix package\n";
  23. //   cout << "\n";
  24.    Tracer et("Seventeenth test of Matrix package");
  25.    Exception::PrintTrace(TRUE);
  26.  
  27.  
  28.    {
  29.       Tracer et1("Stage 1");
  30.       int i, j;
  31.       BandMatrix B(8,3,1);
  32.       for (i=1; i<=8; i++) for (j=-3; j<=1; j++)
  33.          { int k = i+j; if (k>0 && k<=8) B(i,k) = i + k/64.0; }
  34.   
  35.       DiagonalMatrix I(8); I = 1; BandMatrix B1; B1 = I;
  36.       UpperTriangularMatrix UT = I; Print(Matrix(B1-UT));
  37.       Print(Matrix(B * B - B * 2 + I - (B - I) * (B - I)));
  38.       Matrix A = B; BandMatrix C; C = B;
  39.       Print(Matrix(B * A - C * 2 + I - (A - I) * (B - I)));
  40.  
  41.       C.ReDimension(8,2,2); C = 0.0; C.Inject(B);
  42.       Matrix X = A.t();
  43.       B1.ReDimension(8,2,2); B1.Inject(X); Print(Matrix(C.t()-B1));
  44.  
  45.       Matrix BI = B.i(); A = A.i()-BI; Clean(A,0.000000001); Print(A);
  46.  
  47.  
  48.       UpperBandMatrix U; U << B; LowerBandMatrix L; L << B;
  49.       DiagonalMatrix D; D << B;
  50.       Print( Matrix(L + (U - D - B)) );
  51.  
  52.  
  53.  
  54.       for (i=1; i<=8; i++)  A.Column(i) << B.Column(i);
  55.       Print(Matrix(A-B));
  56.    }
  57.    {
  58.       Tracer et1("Stage 2");
  59.       BandMatrix A(7,2,2);
  60.       int i,j;
  61.       for (i=1; i<=7; i++) for (j=1; j<=7; j++)
  62.       {
  63.          int k=i-j; if (k<0) k = -k;
  64.          if (k==0) A(i,j)=6;
  65.          else if (k==1) A(i,j) = -4;
  66.          else if (k==2) A(i,j) = 1;
  67.          A(1,1) = A(7,7) = 5;
  68.       }
  69.       DiagonalMatrix D(7); D = 0.0; A = A - D;
  70.       BandLUMatrix B(A); Matrix M = A;
  71.       ColumnVector V(3);
  72.       V(1) = LogDeterminant(B).Value();
  73.       V(2) = LogDeterminant(A).Value();
  74.       V(3) = LogDeterminant(M).Value();
  75.       V = V / 64 - 1; Clean(V,0.000000001); Print(V);
  76.       ColumnVector X(7);
  77. #ifdef ATandT
  78.       Real a[7]; a[0]=1; a[1]=2; a[2]=3; a[3]=4; a[4]=5; a[5]=6; a[6]=7;
  79. #else
  80.       Real a[] = {1,2,3,4,5,6,7};
  81. #endif
  82.       X << a;
  83.       M = (M.i()*X).t() - (B.i()*X).t() * 2.0 + (A.i()*X).t();
  84.       Clean(M,0.000000001); Print(M);
  85.  
  86.  
  87.       BandMatrix P(80,2,5); ColumnVector CX(80);
  88.       for (i=1; i<=80; i++) for (j=1; j<=80; j++)
  89.       { int d = i-j; if (d<=2 && d>=-5) P(i,j) = i + j/100.0; }
  90.       for (i=1; i<=80; i++)  CX(i) = i*100.0;
  91.       Matrix MP = P;
  92.       ColumnVector V1 = P.i() * CX; ColumnVector V2 = MP.i() * CX;
  93.       V = V1 - V2; Clean(V,0.000000001); Print(V);
  94.  
  95.       V1 = P * V1; V2 = MP * V2; V = V1 - V2; Clean(V,0.000000001); Print(V);
  96.       RowVector XX(1);
  97.       XX = LogDeterminant(P).Value() / LogDeterminant(MP).Value() - 1.0;
  98.       Clean(XX,0.000000001); Print(XX);
  99.  
  100.       LowerBandMatrix LP(80,5);
  101.       for (i=1; i<=80; i++) for (j=1; j<=80; j++)
  102.       { int d = i-j; if (d<=5 && d>=0) LP(i,j) = i + j/100.0; }
  103.       MP = LP;
  104.       XX.ReDimension(4);
  105.       XX(1) = LogDeterminant(LP).Value();
  106.       XX(2) = LogDeterminant(MP).Value();
  107.       V1 = LP.i() * CX; V2 = MP.i() * CX;
  108.       V = V1 - V2; Clean(V,0.000000001); Print(V);
  109.  
  110.       UpperBandMatrix UP(80,4);
  111.       for (i=1; i<=80; i++) for (j=1; j<=80; j++)
  112.       { int d = i-j; if (d<=0 && d>=-4) UP(i,j) = i + j/100.0; }
  113.       MP = UP;
  114.       XX(3) = LogDeterminant(UP).Value();
  115.       XX(4) = LogDeterminant(MP).Value();
  116.       V1 = UP.i() * CX; V2 = MP.i() * CX;
  117.       V = V1 - V2; Clean(V,0.000000001); Print(V);
  118.       XX = XX / SumAbsoluteValue(XX) - .25; Clean(XX,0.000000001); Print(XX);
  119.    }
  120.    {
  121.       Tracer et1("Stage 3");
  122.       SymmetricBandMatrix SA(8,5);
  123.       int i,j;
  124.       for (i=1; i<=8; i++) for (j=1; j<=8; j++)
  125.      if (i-j<=5 && 0<=i-j) SA(i,j) =i + j/128.0;
  126.       DiagonalMatrix D(8); D = 10; SA = SA + D;
  127.  
  128.       Matrix MA1(8,8); Matrix MA2(8,8);
  129.       for (i=1; i<=8; i++)
  130.          { MA1.Column(i) << SA.Column(i); MA2.Row(i) << SA.Row(i); }
  131.       Print(Matrix(MA1-MA2));
  132.  
  133.       D = 10; SA = SA.t() + D; MA1 = MA1 + D;
  134.       Print(Matrix(MA1-SA));
  135.  
  136.       UpperBandMatrix UB(8,3); LowerBandMatrix LB(8,4);
  137.       D << SA; UB << SA; LB << SA;
  138.       SA = SA * 5.0; D = D * 5.0; LB = LB * 5.0; UB = UB * 5.0;
  139.       BandMatrix B = LB - D + UB - SA; Print(Matrix(B));
  140.  
  141.       SymmetricBandMatrix A(7,2); A = 100.0;
  142.       for (i=1; i<=7; i++) for (j=1; j<=7; j++)
  143.       {
  144.      int k=i-j;
  145.      if (k==0) A(i,j)=6;
  146.      else if (k==1) A(i,j) = -4;
  147.      else if (k==2) A(i,j) = 1;
  148.      A(1,1) = A(7,7) = 5;
  149.       }
  150.       BandLUMatrix C(A); Matrix M = A;
  151.       ColumnVector X(7);
  152.       X(1) = LogDeterminant(C).Value() - 64;
  153.       X(2) = LogDeterminant(A).Value() - 64;
  154.       X(3) = LogDeterminant(M).Value() - 64;
  155.       X(4) = SumSquare(M) - SumSquare(A);
  156.       X(5) = SumAbsoluteValue(M) - SumAbsoluteValue(A);
  157.       X(6) = MaximumAbsoluteValue(M) - MaximumAbsoluteValue(A);
  158.       X(7) = Trace(M) - Trace(A);
  159.       Clean(X,0.000000001); Print(X);
  160.  
  161. #ifdef ATandT
  162.       Real a[7]; a[0]=1; a[1]=2; a[2]=3; a[3]=4; a[4]=5; a[5]=6; a[6]=7;
  163. #else
  164.       Real a[] = {1,2,3,4,5,6,7};
  165. #endif
  166.       X << a;
  167.       X = M.i()*X - C.i()*X * 2 + A.i()*X;
  168.       Clean(X,0.000000001); Print(X);
  169.  
  170.  
  171.       LB << A; UB << A; D << A;
  172.       BandMatrix XA = LB + (UB - D);
  173.       Print(Matrix(XA - A));
  174.  
  175.       for (i=1; i<=7; i++) for (j=1; j<=7; j++)
  176.       {
  177.      int k=i-j;
  178.      if (k==0) A(i,j)=6;
  179.      else if (k==1) A(i,j) = -4;
  180.      else if (k==2) A(i,j) = 1;
  181.      A(1,1) = A(7,7) = 5;
  182.       }
  183.       D = 1;
  184.  
  185.       M = LB.i() * LB - D; Clean(M,0.000000001); Print(M);
  186.       M = UB.i() * UB - D; Clean(M,0.000000001); Print(M);
  187.       M = XA.i() * XA - D; Clean(M,0.000000001); Print(M);
  188. #ifdef GXX
  189.       Matrix MUB = UB; Matrix MLB = LB;
  190.       M = LB.i() * UB - LB.i() * MUB; Clean(M,0.000000001); Print(M);
  191.       M = UB.i() * LB - UB.i() * MLB; Clean(M,0.000000001); Print(M);
  192. #else
  193.       M = LB.i() * UB - LB.i() * Matrix(UB); Clean(M,0.000000001); Print(M);
  194.       M = UB.i() * LB - UB.i() * Matrix(LB); Clean(M,0.000000001); Print(M);
  195. #endif
  196.    }
  197.  
  198.  
  199. //   cout << "\nEnd of Seventeenth test\n";
  200. }
  201.