tmtm.cpp
Go to the documentation of this file.
1 
6 
7 
8 #define WANT_STREAM
9 
10 #define WANT_MATH
11 
12 #include "newmat.h"
13 #include "newmatio.h"
14 
15 #include "tmt.h"
16 
17 #ifdef use_namespace
18 using namespace NEWMAT;
19 #endif
20 
21 
22 
23 // test Kronecker Product
24 
25 
26 void trymatm()
27 {
28  Tracer et("Twenty second test of Matrix package");
30 
31  {
32  Tracer et1("Stage 1");
33 
34 
35  Matrix A(2,3);
36  A << 3 << 5 << 2
37  << 4 << 1 << 6;
38 
39  Matrix B(4,3);
40  B << 7 << 2 << 9
41  << 1 << 3 << 6
42  << 4 << 10 << 5
43  << 11 << 8 << 12;
44 
45  Matrix C(8, 9);
46 
47  C.Row(1) << 21 << 6 << 27 << 35 << 10 << 45 << 14 << 4 << 18;
48  C.Row(2) << 3 << 9 << 18 << 5 << 15 << 30 << 2 << 6 << 12;
49  C.Row(3) << 12 << 30 << 15 << 20 << 50 << 25 << 8 << 20 << 10;
50  C.Row(4) << 33 << 24 << 36 << 55 << 40 << 60 << 22 << 16 << 24;
51 
52  C.Row(5) << 28 << 8 << 36 << 7 << 2 << 9 << 42 << 12 << 54;
53  C.Row(6) << 4 << 12 << 24 << 1 << 3 << 6 << 6 << 18 << 36;
54  C.Row(7) << 16 << 40 << 20 << 4 << 10 << 5 << 24 << 60 << 30;
55  C.Row(8) << 44 << 32 << 48 << 11 << 8 << 12 << 66 << 48 << 72;
56 
57  Matrix AB = KP(A,B) - C; Print(AB);
58 
59  IdentityMatrix I1(10); IdentityMatrix I2(15); I2 *= 2;
60  DiagonalMatrix D = KP(I1, I2) - IdentityMatrix(150) * 2;
61  Print(D);
62  }
63 
64  {
65  Tracer et1("Stage 2");
66 
68  A << 3 << 8 << 5
69  << 7 << 2
70  << 4;
72  B << 4 << 1 << 7 << 2
73  << 3 << 9 << 8
74  << 1 << 5
75  << 6;
76 
78 
79  C.Row(1) <<12<< 3<<21<< 6 <<32<< 8<<56<<16 <<20<< 5<<35<<10;
80  C.Row(2) << 9<<27<<24 << 0<<24<<72<<64 << 0<<15<<45<<40;
81  C.Row(3) << 3<<15 << 0<< 0<< 8<<40 << 0<< 0<< 5<<25;
82  C.Row(4) <<18 << 0<< 0<< 0<<48 << 0<< 0<< 0<<30;
83 
84  C.Row(5) <<28<< 7<<49<<14 << 8<< 2<<14<< 4;
85  C.Row(6) <<21<<63<<56 << 0<< 6<<18<<16;
86  C.Row(7) << 7<<35 << 0<< 0<< 2<<10;
87  C.Row(8) <<42 << 0<< 0<< 0<<12;
88 
89  C.Row(9) <<16<< 4<<28<< 8;
90  C.Row(10) <<12<<36<<32;
91  C.Row(11) << 4<<20;
92  C.Row(12) <<24;
93 
94 
95  UpperTriangularMatrix AB = KP(A,B) - C; Print(AB);
96 
97  LowerTriangularMatrix BT = B.t(); Matrix N(12,12);
98 
99  N.Row(1) <<12 << 0<< 0<< 0 <<32<< 0<< 0<< 0 <<20<< 0<< 0<< 0;
100  N.Row(2) << 3 << 9<< 0<< 0 << 8<<24<< 0<< 0 << 5<<15<< 0<< 0;
101  N.Row(3) <<21 <<27<< 3<< 0 <<56<<72<< 8<< 0 <<35<<45<< 5<< 0;
102  N.Row(4) << 6 <<24<<15<<18 <<16<<64<<40<<48 <<10<<40<<25<<30;
103 
104  N.Row(5) << 0 << 0<< 0<< 0 <<28<< 0<< 0<< 0 << 8<< 0<< 0<< 0;
105  N.Row(6) << 0 << 0<< 0<< 0 << 7<<21<< 0<< 0 << 2<< 6<< 0<< 0;
106  N.Row(7) << 0 << 0<< 0<< 0 <<49<<63<< 7<< 0 <<14<<18<< 2<< 0;
107  N.Row(8) << 0 << 0<< 0<< 0 <<14<<56<<35<<42 << 4<<16<<10<<12;
108 
109  N.Row(9) << 0 << 0<< 0<< 0 << 0<< 0<< 0<< 0 <<16<< 0<< 0<< 0;
110  N.Row(10)<< 0 << 0<< 0<< 0 << 0<< 0<< 0<< 0 << 4<<12<< 0<< 0;
111  N.Row(11)<< 0 << 0<< 0<< 0 << 0<< 0<< 0<< 0 <<28<<36<< 4<< 0;
112  N.Row(12)<< 0 << 0<< 0<< 0 << 0<< 0<< 0<< 0 << 8<<32<<20<<24;
113 
114  Matrix N1 = KP(A, BT); N1 -= N; Print(N1);
115  AB << KP(A, BT); AB << (AB - N); Print(AB);
116  BT << KP(A, BT); BT << (BT - N); Print(BT);
117 
118  LowerTriangularMatrix AT = A.t();
119  N1 = KP(AT, B); N1 -= N.t(); Print(N1);
120  AB << KP(AT, B); AB << (AB - N.t()); Print(AB);
121  BT << KP(AT, B); BT << (BT - N.t()); Print(BT);
122  }
123 
124  {
125  Tracer et1("Stage 3");
126 
127  BandMatrix BMA(6,2,3);
128  BMA.Row(1) << 5.25 << 4.75 << 2.25 << 1.75;
129  BMA.Row(2) << 1.25 << 9.75 << 4.50 << 0.25 << 1.50;
130  BMA.Row(3) << 7.75 << 1.50 << 3.00 << 4.25 << 0.50 << 5.50;
131  BMA.Row(4) << 2.75 << 9.00 << 8.00 << 3.25 << 3.50;
132  BMA.Row(5) << 8.75 << 6.25 << 5.00 << 5.75;
133  BMA.Row(6) << 3.75 << 6.75 << 6.00;
134 
135  Matrix A = BMA;
136 
137  BandMatrix BMB(4,2,1);
138  BMB.Row(1) << 4.5 << 9.5;
139  BMB.Row(2) << 1.5 << 6.0 << 2.0;
140  BMB.Row(3) << 0.5 << 2.5 << 8.5 << 7.5;
141  BMB.Row(4) << 3.0 << 4.0 << 6.5;
142 
143  SquareMatrix B = BMB;
144 
145  BandMatrix BMC = KP(BMA, BMB);
146  BandMatrix BMC1 = KP(BMA, B);
147  Matrix C2 = KP(A, BMB);
148  Matrix C = KP(A, B);
149 
150  Matrix M = C - BMC; Print(M);
151  M = C - BMC1; Print(M);
152  M = C - C2; Print(M);
153 
154  RowVector X(4);
155  X(1) = BMC.BandWidth().Lower() - 10;
156  X(2) = BMC.BandWidth().Upper() - 13;
157  X(3) = BMC1.BandWidth().Lower() - 11;
158  X(4) = BMC1.BandWidth().Upper() - 15;
159  Print(X);
160 
161  UpperTriangularMatrix UT; UT << KP(BMA, BMB);
162  UpperTriangularMatrix UT1; UT1 << (C - UT); Print(UT1);
163  LowerTriangularMatrix LT; LT << KP(BMA, BMB);
164  LowerTriangularMatrix LT1; LT1 << (C - LT); Print(LT1);
165  }
166 
167  {
168  Tracer et1("Stage 4");
169 
170  SymmetricMatrix SM1(4);
171  SM1.Row(1) << 2;
172  SM1.Row(2) << 4 << 5;
173  SM1.Row(3) << 9 << 2 << 1;
174  SM1.Row(4) << 3 << 6 << 8 << 2;
175 
176  SymmetricMatrix SM2(3);
177  SM2.Row(1) << 3;
178  SM2.Row(2) << -7 << -6;
179  SM2.Row(3) << 4 << -2 << -1;
180 
181  SymmetricMatrix SM = KP(SM1, SM2);
182  Matrix M1 = SM1; Matrix M2 = SM2;
183  Matrix M = KP(SM1, SM2); M -= SM; Print(M);
184  M = KP(SM1, SM2) - SM; Print(M);
185  M = KP(M1, SM2) - SM; Print(M);
186  M = KP(SM1, M2) - SM; Print(M);
187  M = KP(M1, M2); M -= SM; Print(M);
188  }
189 
190  {
191  Tracer et1("Stage 5");
192 
193  Matrix A(2,3);
194  A << 3 << 5 << 2
195  << 4 << 1 << 6;
196 
197  Matrix B(3,4);
198  B << 7 << 2 << 9 << 11
199  << 1 << 3 << 6 << 8
200  << 4 << 10 << 5 << 12;
201 
202  RowVector C(2); C << 3 << 7;
203  ColumnVector D(4); D << 0 << 5 << 13 << 11;
204 
205  Matrix M = KP(C * A, B * D) - KP(C, B) * KP(A, D); Print(M);
206  }
207 
208  {
209  Tracer et1("Stage 6");
210 
211  RowVector A(3), B(5), C(15);
212  A << 5 << 2 << 4;
213  B << 3 << 2 << 0 << 1 << 6;
214  C << 15 << 10 << 0 << 5 << 30
215  << 6 << 4 << 0 << 2 << 12
216  << 12 << 8 << 0 << 4 << 24;
217  Matrix N = KP(A, B) - C; Print(N);
218  N = KP(A.t(), B.t()) - C.t(); Print(N);
219  N = KP(A.AsDiagonal(), B.AsDiagonal()) - C.AsDiagonal(); Print(N);
220  }
221 
222  {
223  Tracer et1("Stage 7");
224  IdentityMatrix I(3);
225  ColumnVector CV(4); CV << 4 << 3 << 1 << 7;
226  Matrix A = KP(I, CV) + 5;
227  Matrix B(3,12);
228  B.Row(1) << 9 << 8 << 6 << 12 << 5 << 5 << 5 << 5 << 5 << 5 << 5 << 5;
229  B.Row(2) << 5 << 5 << 5 << 5 << 9 << 8 << 6 << 12 << 5 << 5 << 5 << 5;
230  B.Row(3) << 5 << 5 << 5 << 5 << 5 << 5 << 5 << 5 << 9 << 8 << 6 << 12;
231  B -= A.t(); Print(B);
232 
233  }
234 
235  {
236  Tracer et1("Stage 8"); // SquareMatrix
237  Matrix A(2,3), B(3,2);
238  A << 2 << 6 << 7
239  << 4 << 3 << 9;
240  B << 1 << 3
241  << 4 << 8
242  << 0 << 6;
243  SquareMatrix AB = A * B;
244  Matrix M = (B.t() * A.t()).t(); M -= AB; Print(M);
245  AB = B * A;
246  M = (A.t() * B.t()).t(); M -= AB; Print(M);
247  AB.ReSize(5,5); AB = 0;
248  AB.SubMatrix(1,2,1,3) = A; AB.SubMatrix(4,5,3,5) = A;
249  AB.SubMatrix(1,3,4,5) = B; AB.SubMatrix(3,5,1,2) = B;
250  SquareMatrix C(5);
251  C.Row(1) << 2 << 6 << 7 << 1 << 3;
252  C.Row(2) << 4 << 3 << 9 << 4 << 8;
253  C.Row(3) << 1 << 3 << 0 << 0 << 6;
254  C.Row(4) << 4 << 8 << 2 << 6 << 7;
255  C.Row(5) << 0 << 6 << 4 << 3 << 9;
256  C -= AB; Print(C);
257  AB = A.SymSubMatrix(1,2);
258  AB = (AB | AB) & (AB | AB);
259  C.ReSize(4);
260  C.Row(1) << 2 << 6 << 2 << 6;
261  C.Row(2) << 4 << 3 << 4 << 3;
262  C.Row(3) << 2 << 6 << 2 << 6;
263  C.Row(4) << 4 << 3 << 4 << 3;
264  M = AB;
265  C -= M; Print(C);
266  C << M; C += -M; Print(C);
267 
268  }
269 
270 
271 }
272 
273 
275 
276 
277 
278 
279 
KPMatrix KP(const BaseMatrix &, const BaseMatrix &)
Definition: newmat6.cpp:281
virtual MatrixBandWidth BandWidth() const
Definition: newmat.h:389
Upper triangular matrix.
Definition: newmat.h:799
Square matrix.
Definition: newmat.h:679
int Upper() const
Definition: newmat.h:216
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
GetSubMatrix SymSubMatrix(int f, int l) const
Definition: newmat.h:2148
Diagonal matrix.
Definition: newmat.h:896
Lower triangular matrix.
Definition: newmat.h:848
void trymatm()
Definition: tmtm.cpp:26
GetSubMatrix Row(int f) const
Definition: newmat.h:2150
GetSubMatrix SubMatrix(int fr, int lr, int fc, int lc) const
Definition: newmat.h:2146
DiagedMatrix AsDiagonal() const
Definition: newmat.h:2143
Row vector.
Definition: newmat.h:953
int Lower() const
Definition: newmat.h:218
void Print(const Matrix &X)
Definition: tmt.cpp:42
Column vector.
Definition: newmat.h:1008
Identity matrix.
Definition: newmat.h:1350
void ReSize(int m)
Definition: newmat.h:695
Symmetric matrix.
Definition: newmat.h:753


kni
Author(s): Martin Günther
autogenerated on Fri Jun 7 2019 22:06:45