tmta.cpp
Go to the documentation of this file.
1 
6 
7 
8 //#define WANT_STREAM
9 
10 
11 #include "include.h"
12 
13 #include "newmatap.h"
14 
15 #include "tmt.h"
16 
17 #ifdef use_namespace
18 using namespace NEWMAT;
19 #endif
20 
21 
22 /**************************** test program ******************************/
23 
24 
25 static void process(const GeneralMatrix& A,
26  const ColumnVector& X1, const ColumnVector& X2)
27 {
28  Matrix B = A;
30  Matrix Y(4,2);
31  Y.Column(1) << L.i() * X1; Y.Column(2) << L.i() * X2;
32  Matrix Z(4,2); Z.Column(1) << X1; Z.Column(2) << X2;
33  Z = B * Y - Z; Clean(Z,0.00000001); Print(Z);
34 }
35 
36 
37 
38 void trymata()
39 {
40 // cout << "\nTenth test of Matrix package\n";
41  Tracer et("Tenth test of Matrix package");
43  int i; int j;
45  for (i=1;i<=8;i++) for (j=i;j<=8;j++) U(i,j)=i+j*j+5;
46  Matrix X(8,6);
47  for (i=1;i<=8;i++) for (j=1;j<=6;j++) X(i,j)=i*j+1.0;
48  Matrix Y = U.i()*X; Matrix MU=U;
49  Y=Y-MU.i()*X; Clean(Y,0.00000001); Print(Y);
50  Y = U.t().i()*X; Y=Y-MU.t().i()*X; Clean(Y,0.00000001); Print(Y);
52  for (i=1;i<=8;i++) for (j=i;j<=8;j++) UX(i,j)=i+j+1;
53  UX(4,4)=0; UX(4,5)=0;
54  UpperTriangularMatrix UY = U.i() * UX;
55  { X=UX; MU=U; Y=UY-MU.i()*X; Clean(Y,0.000000001); Print(Y); }
56  LowerTriangularMatrix LY = U.t().i() * UX.t();
57  { Y=LY-MU.i().t()*X.t(); Clean(Y,0.000000001); Print(Y); }
58  DiagonalMatrix D(8); for (i=1;i<=8;i++) D(i,i)=i+1;
59  { X=D.i()*MU; }
60  { UY=U; UY=D.i()*UY; Y=UY-X; Clean(Y,0.00000001); Print(Y); }
61  { UY=D.i()*U; Y=UY-X; Clean(Y,0.00000001); Print(Y); }
62 // X=MU.t();
63 // LY=D.i()*U.t(); Y=D*LY-X; Clean(Y,0.00000001); Print(Y);
64 // LowerTriangularMatrix L=U.t();
65 // LY=D.i()*L; Y=D*LY-X; Clean(Y,0.00000001); Print(Y);
66  U.ReSize(8);
67  for (i=1;i<=8;i++) for (j=i;j<=8;j++) U(i,j)=i+j*j+5;
68  MU = U;
69  MU = U.i() - MU.i(); Clean(MU,0.00000001); Print(MU);
70  MU = U.t().i() - U.i().t(); Clean(MU,0.00000001); Print(MU);
71 
72  // test LINEQ
73  {
74  ColumnVector X1(4), X2(4);
75  X1(1)=1; X1(2)=2; X1(3)=3; X1(4)=4;
76  X2(1)=1; X2(2)=10; X2(3)=100; X2(4)=1000;
77 
78 
79  Matrix A(4,4);
80  A(1,1)=1; A(1,2)=3; A(1,3)=0; A(1,4)=0;
81  A(2,1)=3; A(2,2)=2; A(2,3)=5; A(2,4)=0;
82  A(3,1)=0; A(3,2)=5; A(3,3)=4; A(3,4)=1;
83  A(4,1)=0; A(4,2)=0; A(4,3)=1; A(4,4)=3;
84  process(A,X1,X2);
85 
86  BandMatrix B(4,1,1); B.Inject(A);
87  process(B,X1,X2);
88 
90  U(1,1)=1; U(1,2)=2; U(1,3)=3; U(1,4)=4;
91  U(2,2)=8; U(2,3)=7; U(2,4)=6;
92  U(3,3)=2; U(3,4)=7;
93  U(4,4)=1;
94  process(U,X1,X2);
95 
96  // check rowwise load
98  U1.Row(1) << 1 << 2 << 3 << 4;
99  U1.Row(2) << 8 << 7 << 6;
100  U1.Row(3) << 2 << 7;
101  U1.Row(4) << 1;
102 
103  U1 -= U;
104 
105  Print(U1);
106 
107  LowerTriangularMatrix L = U.t();
108  process(L,X1,X2);
109  }
110 
111  // test inversion of poorly conditioned matrix
112  // a user complained this didn't work under OS9
113  {
114  Matrix M(4,4);
115 
116  M << 8.613057e+00 << 8.693985e+00 << -2.322050e-01 << 0.000000e+00
117  << 8.693985e+00 << 8.793868e+00 << -2.346310e-01 << 0.000000e+00
118  << -2.322050e-01 << -2.346310e-01 << 6.264000e-03 << 0.000000e+00
119  << 0.000000e+00 << 0.000000e+00 << 0.000000e+00 << 3.282806e+03 ;
120  Matrix MI = M.i();
121  DiagonalMatrix I(4); I = 1;
122  Matrix Diff = MI * M - I;
123  Clean(Diff,0.00000001); Print(Diff);
124  // Alternatively do Cholesky
125  SymmetricMatrix SM; SM << M;
126  LowerTriangularMatrix LT = Cholesky(SM).i();
127  MI = LT.t() * LT; Diff = MI * M - I;
128  Clean(Diff,0.00000001); Print(Diff);
129  }
130 
131 // cout << "\nEnd of tenth test\n";
132 }
133 
134 
void Inject(const GeneralMatrix &GM)
Definition: newmat.h:526
Upper triangular matrix.
Definition: newmat.h:799
void Clean(Matrix &A, Real c)
Definition: tmt.cpp:128
GetSubMatrix Column(int f) const
Definition: newmat.h:2152
Band matrix.
Definition: newmat.h:1096
TransposedMatrix t() const
Definition: newmat6.cpp:320
static void process(const GeneralMatrix &A, const ColumnVector &X1, const ColumnVector &X2)
Definition: tmta.cpp:25
static void PrintTrace()
Definition: myexcept.cpp:109
The usual rectangular matrix.
Definition: newmat.h:625
InvertedMatrix i() const
Definition: newmat6.cpp:329
ReturnMatrix Cholesky(const SymmetricMatrix &S)
Definition: cholesky.cpp:36
Diagonal matrix.
Definition: newmat.h:896
Lower triangular matrix.
Definition: newmat.h:848
GetSubMatrix Row(int f) const
Definition: newmat.h:2150
void Print(const Matrix &X)
Definition: tmt.cpp:42
Column vector.
Definition: newmat.h:1008
The classes for matrices that can contain data are derived from this.
Definition: newmat.h:447
void trymata()
Definition: tmta.cpp:38
void ReSize(int m)
Definition: newmat.h:833
Symmetric matrix.
Definition: newmat.h:753


kni
Author(s): Martin Günther
autogenerated on Fri Jan 3 2020 04:01:17