tmt9.cpp
Go to the documentation of this file.
1 
6 
7 //#define WANT_STREAM
8 
9 #include "include.h"
10 #include "newmatap.h"
11 
12 #include "tmt.h"
13 
14 #ifdef use_namespace
15 using namespace NEWMAT;
16 #endif
17 
18 
19 /**************************** test program ******************************/
20 
21 
22 void trymat9()
23 {
24  Tracer et("Ninth test of Matrix package");
26 
27 
28  int i; int j;
29  Matrix A(7,7); Matrix X(7,3);
30  for (i=1;i<=7;i++) for (j=1;j<=7;j++) A(i,j)=i*i+j+((i==j) ? 1 : 0);
31  for (i=1;i<=7;i++) for (j=1;j<=3;j++) X(i,j)=i-j;
32  Matrix B = A.i(); DiagonalMatrix D(7); D=1.0;
33  {
34  Tracer et1("Stage 1");
35  Matrix Q = B*A-D; Clean(Q, 0.000000001); Print(Q);
36  Q=A; Q = Q.i() * X; Q = A*Q - X; Clean(Q, 0.000000001); Print(Q);
37  Q=X; Q = A.i() * Q; Q = A*Q - X; Clean(Q, 0.000000001); Print(Q);
38  }
39  for (i=1;i<=7;i++) D(i,i)=i*i+1;
40  DiagonalMatrix E(3); for (i=1;i<=3;i++) E(i,i)=i+23;
41  {
42  Tracer et1("Stage 2");
43  Matrix DXE = D.i() * X * E;
44  DXE = E.i() * DXE.t() * D - X.t(); Clean(DXE, 0.00000001); Print(DXE);
45  E=D; for (i=1;i<=7;i++) E(i,i)=i*3+1;
46  }
47  DiagonalMatrix F=D;
48  {
49  Tracer et1("Stage 3");
50  F=E.i()*F; F=F*E-D; Clean(F,0.00000001); Print(F);
51  F=E.i()*D; F=F*E-D; Clean(F,0.00000001); Print(F);
52  }
53  {
54  Tracer et1("Stage 4");
55  F=E; F=F.i()*D; F=F*E-D; Clean(F,0.00000001); Print(F);
56  }
57  {
58  Tracer et1("Stage 5");
59  // testing equal
60  ColumnVector A(18), B(18);
61  Matrix X(3,3);
62  X << 3 << 5 << 7 << 5 << 8 << 2 << 7 << 2 << 9;
63  SymmetricMatrix S; S << X;
64  B(1) = S == X; A(1) = true;
65  B(2) = S == (X+1); A(2) = false;
66  B(3) = (S+2) == (X+2); A(3) = true;
67  Matrix Y = X;
68  B(4) = X == Y; A(4) = true;
69  B(5) = (X*2) == (Y*2); A(5) = true;
70  Y(3,3) = 10;
71  B(6) = X == Y; A(6) = false;
72  B(7) = (X*2) == (Y*2); A(7) = false;
73  B(8) = S == Y; A(8) = false;
74  B(9) = S == S; A(9) = true;
75  Matrix Z = X.SubMatrix(1,2,2,3);
76  B(10) = X == Z; A(10) = false;
77  GenericMatrix GS = S;
78  GenericMatrix GX = X;
79  GenericMatrix GY = Y;
80  B(11) = GS == GX; A(11) = true;
81  B(12) = GS == GY; A(12) = false;
82  CroutMatrix CS = S;
83  CroutMatrix CX = X;
84  CroutMatrix CY = Y;
85  B(13) = CS == CX; A(13) = true;
86  B(14) = CS == CY; A(14) = false;
87  B(15) = X == CX; A(15) = false;
88  B(16) = X == A; A(16) = false;
89  B(17) = X == (X | X); A(17) = false;
90  B(18) = CX == X; A(18) = false;
91  A = A - B; Print(A);
92  }
93  {
94  Tracer et1("Stage 6");
95  // testing equal
96  ColumnVector A(22), B(22);
97  BandMatrix X(6,2,1);
98  X(1,1)=23; X(1,2)=21;
99  X(2,1)=12; X(2,2)=17; X(2,3)=45;
100  X(3,1)=35; X(3,2)=19; X(3,3)=24; X(3,4)=29;
101  X(4,2)=17; X(4,3)=11; X(4,4)=19; X(4,5)=35;
102  X(5,3)=10; X(5,4)=44; X(5,5)=23; X(5,6)=31;
103  X(6,4)=49; X(6,5)=41; X(6,6)=17;
104  SymmetricBandMatrix S1(6,2); S1.Inject(X);
105  BandMatrix U(6,2,3); U = 0.0; U.Inject(X);
106  B(1) = U == X; A(1) = true;
107  B(2) = U == (X*3); A(2) = false;
108  B(3) = (U*5) == (X*5); A(3) = true;
109  Matrix Y = X;
110  B(4) = X == Y; A(4) = true;
111  B(5) = (X*2) == (Y*2); A(5) = true;
112  Y(6,6) = 10;
113  B(6) = X == Y; A(6) = false;
114  B(7) = (X*2) == (Y*2); A(7) = false;
115  B(8) = U == Y; A(8) = false;
116  B(9) = U == U; A(9) = true;
117  Matrix Z = X.SubMatrix(1,2,2,3);
118  B(10) = X == Z; A(10) = false;
119  GenericMatrix GU = U;
120  GenericMatrix GX = X;
121  GenericMatrix GY = Y;
122  B(11) = GU == GX; A(11) = true;
123  B(12) = GU == GY; A(12) = false;
124  X = X + X.t(); U = U + U.t();
125  SymmetricBandMatrix S(6,2); S.Inject(X);
126  Matrix D = S-X; Print(D);
127  BandLUMatrix BS = S;
128  BandLUMatrix BX = X;
129  BandLUMatrix BU = U;
130  CroutMatrix CX = X;
131  B(13) = BS == BX; A(13) = true;
132  B(14) = BX == BU; A(14) = false;
133  B(15) = X == BX; A(15) = false;
134  B(16) = X != BX; A(16) = true;
135  B(17) = BX != BS; A(17) = false;
136  B(18) = (2*X) != (X*2);A(18) = false;
137  B(19) = (X*2) != (X+2);A(19) = true;
138  B(20) = BX == CX; A(20) = false;
139  B(21) = CX == BX; A(21) = false;
140  B(22) = BX == X; A(22) = false;
141  A = A - B; Print(A);
142  DiagonalMatrix I(6); I=1.0;
143  D = BS.i() * X - I; Clean(D,0.00000001); Print(D);
144  D = BX.i() * X - I; Clean(D,0.00000001); Print(D);
145  D = BU.i() * X - I; Clean(D,0.00000001); Print(D);
146 
147  // test row wise load
148  SymmetricBandMatrix X1(6,2);
149  X1.Row(1) << 23;
150  X1.Row(2) << 12 << 17;
151  X1.Row(3) << 35 << 19 << 24;
152  X1.Row(4) << 17 << 11 << 19;
153  X1.Row(5) << 10 << 44 << 23;
154  X1.Row(6) << 49 << 41 << 17;
155  Matrix M = X1 - S1; Print(M);
156 
157  // check out submatrix
158  SymmetricBandMatrix X2(20,3); X2 = 0.0;
159  X2.SubMatrix(2,7,2,7) = X1; X2.SymSubMatrix(11,16) = 2 * X1;
160  Matrix MX1 = X1;
161  Matrix MX2(20,20); MX2 = 0;
162  MX2.SymSubMatrix(2,7) = MX1; MX2.SubMatrix(11,16,11,16) = MX1 * 2;
163  MX2 -= X2; Print(MX2);
164 
165  BandMatrix X4(20,3,3); X4 = 0.0;
166  X4.SubMatrix(2,7,3,8) = X1; X4.SubMatrix(11,16,10,15) = 2 * X1;
167  MX1 = X1;
168  Matrix MX4(20,20); MX4 = 0;
169  MX4.SubMatrix(2,7,3,8) = MX1; MX4.SubMatrix(11,16,10,15) = MX1 * 2;
170  MX4 -= X4; Print(MX4);
171 
172  MX1 = X1.i() * X1 - IdentityMatrix(6);
173  Clean(MX1,0.00000001); Print(MX1);
174 
175  }
176 
177  {
178  Tracer et1("Stage 7");
179  // testing equal
180  ColumnVector A(12), B(12);
181  BandMatrix X(6,2,1);
182  X(1,1)=23; X(1,2)=21;
183  X(2,1)=12; X(2,2)=17; X(2,3)=45;
184  X(3,1)=35; X(3,2)=19; X(3,3)=24; X(3,4)=29;
185  X(4,2)=17; X(4,3)=11; X(4,4)=19; X(4,5)=35;
186  X(5,3)=10; X(5,4)=44; X(5,5)=23; X(5,6)=31;
187  X(6,4)=49; X(6,5)=41; X(6,6)=17;
188  Matrix Y = X;
189  LinearEquationSolver LX = X;
190  LinearEquationSolver LY = Y;
191  CroutMatrix CX = X;
192  CroutMatrix CY = Y;
193  BandLUMatrix BX = X;
194  B(1) = LX == CX; A(1) = false;
195  B(2) = LY == CY; A(2) = true;
196  B(3) = X == Y; A(3) = true;
197  B(4) = BX == LX; A(4) = true;
198  B(5) = CX == CY; A(5) = true;
199  B(6) = LX == LY; A(6) = false;
200  B(7) = BX == BX; A(7) = true;
201  B(8) = CX == CX; A(8) = true;
202  B(9) = LX == LX; A(9) = true;
203  B(10) = LY == LY; A(10) = true;
204  CroutMatrix CX1 = X.SubMatrix(1,4,1,4);
205  B(11) = CX == CX1; A(11) = false;
206  BandLUMatrix BX1 = X.SymSubMatrix(1,4); // error with SubMatrix
207  B(12) = BX == BX1; A(12) = false;
208  A = A - B; Print(A);
209  DiagonalMatrix I(6); I=1.0; Matrix D;
210  D = LX.i() * X - I; Clean(D,0.00000001); Print(D);
211  D = LY.i() * X - I; Clean(D,0.00000001); Print(D);
212  I.ReSize(4); I = 1;
213  D = CX1.i() * X.SymSubMatrix(1,4) - I; Clean(D,0.00000001); Print(D);
214  D = BX1.i() * X.SubMatrix(1,4,1,4) - I; Clean(D,0.00000001); Print(D);
215  }
216 
217  {
218  Tracer et1("Stage 8");
219  // test copying CroutMatrix and BandLUMatrix - see also tmtd.cpp
220  MultWithCarry MWC;
221  SymmetricBandMatrix SBM(50, 10);
222  for (int i = 1; i <= 50; ++i) for (int j = 1; j <= i; ++j)
223  if (i - j <= 10) SBM(i, j) = MWC.Next();
224  CroutMatrix CM = SBM; BandLUMatrix BM = SBM;
225  CroutMatrix CM1 = CM; BandLUMatrix BM1; BM1 = BM;
226  CM1.release(); BM1.release();
227  CroutMatrix CM2; CM2 = CM1; BandLUMatrix BM2 = BM1;
228  Matrix X = SBM.i(); Matrix Y = CM2.i() - X; Matrix Z = BM2.i() - X;
229  Clean(Y,0.00000001); Print(Y); Clean(Z,0.00000001); Print(Z);
230  X *= SBM; X -= IdentityMatrix(50); Clean(X,0.00000001); Print(X);
231  LogAndSign x = log_determinant(SBM);
232  LogAndSign y = log_determinant(CM2);
233  LogAndSign z = log_determinant(BM2);
234  RowVector D(4);
235  D(1) = y.value() - x.value();
236  D(2) = z.value() - x.value();
237  D(3) = y.sign() - x.sign();
238  D(4) = z.sign() - x.sign();
239  Clean(D,0.00000001); Print(D);
240  }
241 
242  {
243  Tracer et1("Stage 9");
244  // do it again odd matrix size
245  MultWithCarry MWC;
246  SymmetricBandMatrix SBM(51, 10);
247  for (int i = 1; i <= 51; ++i) for (int j = 1; j <= i; ++j)
248  if (i - j <= 10) SBM(i, j) = MWC.Next();
249  CroutMatrix CM = SBM; BandLUMatrix BM = SBM;
250  CroutMatrix CM1 = CM; BandLUMatrix BM1; BM1 = BM;
251  CM1.release(); BM1.release();
252  CroutMatrix CM2; CM2 = CM1; BandLUMatrix BM2 = BM1;
253  Matrix X = SBM.i(); Matrix Y = CM2.i() - X; Matrix Z = BM2.i() - X;
254  Clean(Y,0.00000001); Print(Y); Clean(Z,0.00000001); Print(Z);
255  X *= SBM; X -= IdentityMatrix(51); Clean(X,0.00000001); Print(X);
256  LogAndSign x = log_determinant(SBM);
257  LogAndSign y = log_determinant(CM2);
258  LogAndSign z = log_determinant(BM2);
259  RowVector D(4);
260  D(1) = y.value() - x.value();
261  D(2) = z.value() - x.value();
262  D(3) = y.sign() - x.sign();
263  D(4) = z.sign() - x.sign();
264  Clean(D,0.00000001); Print(D);
265  }
266 
267 
268 // cout << "\nEnd of ninth test\n";
269 }
270 
271 
LogAndSign log_determinant(const BaseMatrix &B)
Definition: newmat.h:2086
void trymat9()
Definition: tmt9.cpp:22
void Inject(const GeneralMatrix &GM)
Definition: newmat.h:526
Real Next()
Definition: tmt.cpp:187
void Clean(Matrix &A, Real c)
Definition: tmt.cpp:128
Real value() const
the value
Definition: newmat8.cpp:641
Band matrix.
Definition: newmat.h:1096
TransposedMatrix t() const
Definition: newmat6.cpp:320
static void PrintTrace()
Definition: myexcept.cpp:109
The usual rectangular matrix.
Definition: newmat.h:625
InvertedMatrix i() const
Definition: newmat6.cpp:329
GetSubMatrix SymSubMatrix(int f, int l) const
Definition: newmat.h:2148
LU decomposition of a band matrix.
Definition: newmat.h:1306
int sign() const
sign of the value
Definition: newmat.h:62
Diagonal matrix.
Definition: newmat.h:896
GetSubMatrix Row(int f) const
Definition: newmat.h:2150
GetSubMatrix SubMatrix(int fr, int lr, int fc, int lc) const
Definition: newmat.h:2146
Symmetric band matrix.
Definition: newmat.h:1245
Row vector.
Definition: newmat.h:953
void Print(const Matrix &X)
Definition: tmt.cpp:42
Column vector.
Definition: newmat.h:1008
A matrix which can be of any GeneralMatrix type.
Definition: newmat.h:1397
void release()
Definition: newmat.h:517
Identity matrix.
Definition: newmat.h:1350
Symmetric matrix.
Definition: newmat.h:753


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