home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #1 / NN_1993_1.iso / spool / comp / sources / misc / 4255 / tmtd.cxx < prev    next >
Encoding:
C/C++ Source or Header  |  1993-01-11  |  3.2 KB  |  122 lines
  1.  
  2. //#define WANT_STREAM
  3.  
  4. #include "include.h"
  5. #include "newmatap.h"
  6.  
  7. void Print(const Matrix& X);
  8. void Print(const UpperTriangularMatrix& X);
  9. void Print(const DiagonalMatrix& X);
  10. void Print(const SymmetricMatrix& X);
  11. void Print(const LowerTriangularMatrix& X);
  12.  
  13. void Clean(Matrix&, Real);
  14.  
  15.  
  16.  
  17.  
  18. void trymatd()
  19. {
  20. //   cout << "\nThirteenth test of Matrix package\n";
  21.    Tracer et("Thirteenth test of Matrix package");
  22.    Exception::PrintTrace(TRUE);
  23.    Matrix X(5,20);
  24.    int i,j;
  25.    for (j=1;j<=20;j++) X(1,j) = j+1;
  26.    for (i=2;i<=5;i++) for (j=1;j<=20; j++) X(i,j) = (long)X(i-1,j) * j % 1001;
  27.    SymmetricMatrix S; S << X * X.t();
  28.    Matrix SM = X * X.t() - S;
  29.    Print(SM);
  30.    LowerTriangularMatrix L = Cholesky(S);
  31.    Matrix Diff = L*L.t()-S; Clean(Diff, 0.000000001);
  32.    Print(Diff);
  33.    LowerTriangularMatrix L1(5);
  34.    HHDecompose(X,L1);
  35.    Diff = L - L1; Clean(Diff,0.000000001); Print(Diff);
  36.  
  37.    ColumnVector C1(5);
  38.    {
  39.       Tracer et1("Stage 1");
  40.       X.ReDimension(5,5);
  41.       for (j=1;j<=5;j++) X(1,j) = j+1;
  42.       for (i=2;i<=5;i++) for (j=1;j<=5; j++)
  43.          X(i,j) = (long)X(i-1,j) * j % 1001;
  44.       for (i=1;i<=5;i++) C1(i) = i*i;
  45.       CroutMatrix A = X;
  46.       ColumnVector C2 = A.i() * C1; C1 = X.i()  * C1;
  47.       X = C1 - C2; Clean(X,0.000000001); Print(X); 
  48.    }
  49.  
  50.    {
  51.       Tracer et1("Stage 2");
  52.       X.ReDimension(7,7);
  53.       for (j=1;j<=7;j++) X(1,j) = j+1;
  54.       for (i=2;i<=7;i++) for (j=1;j<=7; j++)
  55.          X(i,j) = (long)X(i-1,j) * j % 1001;
  56.       C1.ReDimension(7);
  57.       for (i=1;i<=7;i++) C1(i) = i*i;
  58.       RowVector R1 = C1.t();
  59.       Diff = R1 * X.i() - ( X.t().i() * R1.t() ).t(); Clean(Diff,0.000000001);
  60.       Print(Diff);
  61.    }
  62.  
  63.    {
  64.       Tracer et1("Stage 3");
  65.       X.ReDimension(5,5);
  66.       for (j=1;j<=5;j++) X(1,j) = j+1;
  67.       for (i=2;i<=5;i++) for (j=1;j<=5; j++)
  68.          X(i,j) = (long)X(i-1,j) * j % 1001;
  69.       C1.ReDimension(5);
  70.       for (i=1;i<=5;i++) C1(i) = i*i;
  71.       CroutMatrix A1 = X*X;
  72.       ColumnVector C2 = A1.i() * C1; C1 = X.i()  * C1; C1 = X.i()  * C1;
  73.       X = C1 - C2; Clean(X,0.000000001); Print(X); 
  74.    }
  75.  
  76.  
  77.    {
  78.       Tracer et1("Stage 4");
  79.       int n = 40;
  80.       SymmetricBandMatrix B(n,2); B = 0.0;
  81.       for (i=1; i<=n; i++)
  82.       {
  83.      B(i,i) = 6;
  84.      if (i<=n-1) B(i,i+1) = -4;
  85.      if (i<=n-2) B(i,i+2) = 1;
  86.       }
  87.       B(1,1) = 5; B(n,n) = 5;
  88.       SymmetricMatrix A = B;
  89.       ColumnVector X(n); X = 0; X(1) = 1;
  90.       ColumnVector Y1 = A.i() * X;
  91.       LowerTriangularMatrix C1 = Cholesky(A);
  92.       ColumnVector Y2 = C1.t().i() * (C1.i() * X) - Y1;
  93.       Clean(Y2, 0.000000001); Print(Y2);
  94.   
  95.       LowerBandMatrix C2 = Cholesky(B);
  96.       Matrix M = C2 - C1; Clean(M, 0.000000001); Print(M);
  97.       ColumnVector Y3 = C2.t().i() * (C2.i() * X) - Y1;
  98.       Clean(Y3, 0.000000001); Print(Y3);
  99.    }
  100.  
  101.    {
  102.       Tracer et1("Stage 5");
  103.       SymmetricMatrix A(4), B(4);
  104.       A << 5
  105.         << 1 << 4
  106.         << 2 << 1 << 6
  107.         << 1 << 0 << 1 << 7;
  108.       B << 8
  109.         << 1 << 5
  110.         << 1 << 0 << 9
  111.         << 2 << 1 << 0 << 6;
  112.       LowerTriangularMatrix AB = Cholesky(A) * Cholesky(B);
  113.       Matrix M = Cholesky(A) * B * Cholesky(A).t() - AB*AB.t();
  114.       Clean(M, 0.000000001); Print(M);
  115.       M = A * Cholesky(B); M = M * M.t() - A * B * A;
  116.       Clean(M, 0.000000001); Print(M);
  117.    }
  118.  
  119.  
  120. //   cout << "\nEnd of Thirteenth test\n";
  121. }
  122.