home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / World_Of_Computer_Software-02-385-Vol-1of3.iso / n / newmat06.zip / TMTD.CXX < prev    next >
C/C++ Source or Header  |  1992-11-26  |  4KB  |  123 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.    LowerTriangularMatrix L = Cholesky(S);
  29.    Matrix Diff = L*L.t()-S; Clean(Diff, 0.000000001);
  30.    Print(Diff);
  31.    LowerTriangularMatrix L1(5);
  32.    HHDecompose(X,L1);
  33.    Diff = L - L1; Clean(Diff,0.000000001); Print(Diff);
  34.  
  35.    ColumnVector C1(5);
  36.    {
  37.       Tracer et1("Stage 1");
  38.       X.ReDimension(5,5);
  39.       for (j=1;j<=5;j++) X(1,j) = j+1;
  40.       for (i=2;i<=5;i++) for (j=1;j<=5; j++)
  41.          X(i,j) = (long)X(i-1,j) * j % 1001;
  42.       for (i=1;i<=5;i++) C1(i) = i*i;
  43.       CroutMatrix A = X;
  44.       ColumnVector C2 = A.i() * C1; C1 = X.i()  * C1;
  45.       X = C1 - C2; Clean(X,0.000000001); Print(X); 
  46.    }
  47.  
  48.    {
  49.       Tracer et1("Stage 2");
  50.       X.ReDimension(7,7);
  51.       for (j=1;j<=7;j++) X(1,j) = j+1;
  52.       for (i=2;i<=7;i++) for (j=1;j<=7; j++)
  53.          X(i,j) = (long)X(i-1,j) * j % 1001;
  54.       C1.ReDimension(7);
  55.       for (i=1;i<=7;i++) C1(i) = i*i;
  56.       RowVector R1 = C1.t();
  57.       Diff = R1 * X.i() - ( X.t().i() * R1.t() ).t(); Clean(Diff,0.000000001);
  58.       Print(Diff);
  59.    }
  60.  
  61.    {
  62.       Tracer et1("Stage 3");
  63.       X.ReDimension(5,5);
  64.       for (j=1;j<=5;j++) X(1,j) = j+1;
  65.       for (i=2;i<=5;i++) for (j=1;j<=5; j++)
  66.          X(i,j) = (long)X(i-1,j) * j % 1001;
  67.       C1.ReDimension(5);
  68.       for (i=1;i<=5;i++) C1(i) = i*i;
  69.       CroutMatrix A1 = X*X;
  70.       ColumnVector C2 = A1.i() * C1; C1 = X.i()  * C1; C1 = X.i()  * C1;
  71.       X = C1 - C2; Clean(X,0.000000001); Print(X); 
  72.    }
  73.  
  74.  
  75.    {
  76.       Tracer et1("Stage 4");
  77.       int n = 40;
  78.       SymmetricBandMatrix B(n,2); B = 0.0;
  79.       for (i=1; i<=n; i++)
  80.       {
  81.      B(i,i) = 6;
  82.      if (i<=n-1) B(i,i+1) = -4;
  83.      if (i<=n-2) B(i,i+2) = 1;
  84.       }
  85.       B(1,1) = 5; B(n,n) = 5;
  86.       SymmetricMatrix A = B;
  87.       ColumnVector X(n); X = 0; X(1) = 1;
  88.       ColumnVector Y1 = A.i() * X;
  89.       LowerTriangularMatrix C1 = Cholesky(A);
  90.       ColumnVector Y2 = C1.t().i() * (C1.i() * X) - Y1;
  91.       Clean(Y2, 0.000000001); Print(Y2);
  92.   
  93.       LowerBandMatrix C2 = Cholesky(B);
  94.       Matrix M = C2 - C1; Clean(M, 0.000000001); Print(M);
  95.       ColumnVector Y3 = C2.t().i() * (C2.i() * X) - Y1;
  96.       Clean(Y3, 0.000000001); Print(Y3);
  97.    }
  98.  
  99.    {
  100.       Tracer et1("Stage 5");
  101.       SymmetricMatrix A(4), B(4);
  102. #ifdef ATandT
  103.       Real a[10], b[10];
  104.       a[0] = 5; a[1] = 1; a[2] = 4; a[3] = 2; a[4] = 1;
  105.       a[5] = 6; a[6] = 1; a[7] = 0; a[8] = 1; a[9] = 7;
  106.       b[0] = 8; b[1] = 1; b[2] = 5; b[3] = 1; b[4] = 0;
  107.       b[5] = 9; b[6] = 2; b[7] = 1; b[8] = 0; b[9] = 6; 
  108. #else
  109.       Real a[] = { 5, 1, 4, 2, 1, 6, 1, 0, 1, 7 };
  110.       Real b[] = { 8, 1, 5, 1, 0, 9, 2, 1, 0, 6 };
  111. #endif
  112.       A << a; B << b;
  113.       LowerTriangularMatrix AB = Cholesky(A) * Cholesky(B);
  114.       Matrix M = Cholesky(A) * B * Cholesky(A).t() - AB*AB.t();
  115.       Clean(M, 0.000000001); Print(M);
  116.       M = A * Cholesky(B); M = M * M.t() - A * B * A;
  117.       Clean(M, 0.000000001); Print(M);
  118.    }
  119.  
  120.  
  121. //   cout << "\nEnd of Thirteenth test\n";
  122. }
  123.