tmt3.cpp
Go to the documentation of this file.
1 
6 
7 
8 //#define WANT_STREAM
9 
10 #include "include.h"
11 
12 #include "newmat.h"
13 
14 #include "tmt.h"
15 
16 #ifdef use_namespace
17 using namespace NEWMAT;
18 #endif
19 
20 
21 /**************************** test program ******************************/
22 
23 void trymat3()
24 {
25  Tracer et("Third test of Matrix package");
27 
28  {
29  Tracer et1("Stage 1");
30  int i,j;
31  SymmetricMatrix S(7);
32  for (i=1;i<=7;i++) for (j=1;j<=i;j++) S(i,j)=i*i+j;
33  S=-S+2.0;
34 
35  DiagonalMatrix D(7);
36  for (i=1;i<=7;i++) D(i,i)=S(i,i);
37 
38  Matrix M4(7,7); { M4=D+(D+4.0); M4=M4-D*2.0; M4=M4-4.0; Print(M4); }
39  SymmetricMatrix S2=D; Matrix M2=S2; { M2=-D+M2; Print(M2); }
40  UpperTriangularMatrix U2=D; { M2=U2; M2=D-M2; Print(M2); }
41  LowerTriangularMatrix L2=D; { M2=L2; M2=D-M2; Print(M2); }
42  M2=D; M2=M2-D; Print(M2);
43  for (i=1;i<=7;i++) for (j=1;j<=i;j++) L2(i,j)=2.0-i*i-j;
44  U2=L2.t(); D=D.t(); S=S.t();
45  M4=(L2-1.0)+(U2+1.0)-D-S; Print(M4);
46  M4=(-L2+1.0)+(-U2-1.0)+D+S; Print(M4);
47  }
48 
49  {
50  Tracer et1("Stage 2");
51  int i,j;
52  DiagonalMatrix D(6);
53  for (i=1;i<=6;i++) D(i,i)=i*3.0+i*i+2.0;
55  for (i=1;i<=7;i++) for (j=1;j<=i;j++) L2(i,j)=2.0-i*i+j;
56  { U2=L2.t(); }
57  DiagonalMatrix D1(7); for (i=1;i<=7;i++) D1(i,i)=(i-2)*(i-4);
58  Matrix M2(6,7);
59  for (i=1;i<=6;i++) for (j=1;j<=7;j++) M2(i,j)=2.0+i*j+i*i+2.0*j*j;
60  Matrix MD=D; SymmetricMatrix MD1(1); MD1=D1;
61  Matrix MX=MD*M2*MD1 - D*(M2*D1); Print(MX);
62  MX=MD*M2*MD1 - (D*M2)*D1; Print(MX);
63  {
64  D.ReSize(7); for (i=1;i<=7;i++) D(i,i)=i*3.0+i*i+2.0;
65  LowerTriangularMatrix LD(1); LD=D;
66  UpperTriangularMatrix UD(1); UD=D;
67  M2=U2; M2=LD*M2*MD1 - D*(U2*D1); Print(M2);
68  M2=U2; M2=UD*M2*MD1 - (D*U2)*D1; Print(M2);
69  M2=L2; M2=LD*M2*MD1 - D*(L2*D1); Print(M2);
70  M2=L2; M2=UD*M2*MD1 - (D*L2)*D1; Print(M2);
71  }
72  }
73 
74  {
75  Tracer et1("Stage 3");
76  // test inverse * scalar
77  DiagonalMatrix D(6);
78  for (int i=1;i<=6;i++) D(i)=i*i;
79  DiagonalMatrix E = D.i() * 4.0;
80  DiagonalMatrix I(6); I = 1.0;
81  E=D*E-I*4.0; Print(E);
82  E = D.i() / 0.25; E=D*E-I*4.0; Print(E);
83  }
84  {
85  Tracer et1("Stage 4");
86  Matrix sigma(3,3); Matrix sigmaI(3,3);
87  sigma = 0; sigma(1,1) = 1.0; sigma(2,2) = 1.0; sigma(3,3) = 1.0;
88  sigmaI = sigma.i();
89  sigmaI -= sigma; Clean(sigmaI, 0.000000001); Print(sigmaI);
90  }
91  {
92  Tracer et1("Stage 5");
93  Matrix X(5,5); DiagonalMatrix DM(5);
94  for (int i=1; i<=5; i++) for (int j=1; j<=5; j++)
95  X(i,j) = (23*i+59*j) % 43;
96  DM << 1 << 8 << -7 << 2 << 3;
97  Matrix Y = X.i() * DM; Y = X * Y - DM;
98  Clean(Y, 0.000000001); Print(Y);
99  }
100  {
101  Tracer et1("Stage 6"); // test reverse function
102  ColumnVector CV(10), RCV(10);
103  CV << 2 << 7 << 1 << 6 << -3 << 1 << 8 << -4 << 0 << 17;
104  RCV << 17 << 0 << -4 << 8 << 1 << -3 << 6 << 1 << 7 << 2;
105  ColumnVector X = CV - RCV.Reverse(); Print(X);
106  RowVector Y = CV.t() - RCV.t().Reverse(); Print(Y);
107  DiagonalMatrix D = CV.AsDiagonal() - RCV.AsDiagonal().Reverse();
108  Print(D);
109  X = CV & CV.Rows(1,9).Reverse();
110  ColumnVector Z(19);
111  Z.Rows(1,10) = RCV.Reverse(); Z.Rows(11,19) = RCV.Rows(2,10);
112  X -= Z; Print(X); Z -= Z.Reverse(); Print(Z);
113  Matrix A(3,3); A << 1 << 2 << 3 << 4 << 5 << 6 << 7 << 8 << 9;
114  Matrix B(3,3); B << 9 << 8 << 7 << 6 << 5 << 4 << 3 << 2 << 1;
115  Matrix Diff = A - B.Reverse(); Print(Diff);
116  Diff = (-A).Reverse() + B; Print(Diff);
118  U << A.Reverse(); Diff = U; U << B; Diff -= U; Print(Diff);
119  U << (-A).Reverse(); Diff = U; U << B; Diff += U; Print(Diff);
120  }
121  {
122  Tracer et1("Stage 7"); // test IsSingular function
123  ColumnVector XX(4);
124  Matrix A(3,3);
125  A = 0;
126  CroutMatrix B1 = A;
127  XX(1) = B1.IsSingular() ? 0 : 1;
128  A << 1 << 3 << 6
129  << 7 << 11 << 13
130  << 2 << 4 << 1;
131  CroutMatrix B2(A);
132  XX(2) = B2.IsSingular() ? 1 : 0;
133  BandMatrix C(3,1,1); C.Inject(A);
134  BandLUMatrix B3(C);
135  XX(3) = B3.IsSingular() ? 1 : 0;
136  C = 0;
137  BandLUMatrix B4(C);
138  XX(4) = B4.IsSingular() ? 0 : 1;
139  Print(XX);
140  }
141  {
142  Tracer et1("Stage 8"); // inverse with vector of 0s
143  Matrix A(3,3); Matrix Z(3,3); ColumnVector X(6);
144  A << 1 << 3 << 6
145  << 7 << 11 << 13
146  << 2 << 4 << 1;
147  Z = 0;
148  Matrix B = (A | Z) & (Z | A); // 6 * 6 matrix
149  X = 0.0;
150  X = B.i() * X;
151  Print(X);
152  // also check inverse with non-zero Y
153  Matrix Y(3,3);
154  Y << 0.0 << 1.0 << 1.0
155  << 5.0 << 0.0 << 5.0
156  << 3.0 << 3.0 << 0.0;
157  Matrix YY = Y & Y; // stack Y matrices
158  YY = B.i() * YY;
159  Matrix Y1 = A.i() * Y;
160  YY -= Y1 & Y1; Clean(YY, 0.000000001); Print(YY);
161  Y1 = A * Y1 - Y; Clean(Y1, 0.000000001); Print(Y1);
162  }
163 
164 
165 }
166 
167 
void Inject(const GeneralMatrix &GM)
Definition: newmat.h:526
void ReSize(int m)
Definition: newmat.h:935
ReversedMatrix Reverse() const
Definition: newmat.h:2140
bool IsSingular() const
Definition: newmat.h:1086
Upper triangular matrix.
Definition: newmat.h:799
void Clean(Matrix &A, Real c)
Definition: tmt.cpp:128
Band matrix.
Definition: newmat.h:1096
TransposedMatrix t() const
Definition: newmat6.cpp:320
bool IsSingular() const
Definition: newmat.h:1337
static void PrintTrace()
Definition: myexcept.cpp:109
The usual rectangular matrix.
Definition: newmat.h:625
InvertedMatrix i() const
Definition: newmat6.cpp:329
LU decomposition of a band matrix.
Definition: newmat.h:1306
Diagonal matrix.
Definition: newmat.h:896
Lower triangular matrix.
Definition: newmat.h:848
DiagedMatrix AsDiagonal() const
Definition: newmat.h:2143
Row vector.
Definition: newmat.h:953
void Print(const Matrix &X)
Definition: tmt.cpp:42
Column vector.
Definition: newmat.h:1008
GetSubMatrix Rows(int f, int l) const
Definition: newmat.h:2151
void trymat3()
Definition: tmt3.cpp:23
Symmetric matrix.
Definition: newmat.h:753


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