tmth.cpp
Go to the documentation of this file.
1 
6 
7 
8 //#define WANT_STREAM
9 
10 #include "include.h"
11 
12 #include "newmatap.h"
13 //#include "newmatio.h"
14 
15 #include "tmt.h"
16 
17 #ifdef use_namespace
18 using namespace NEWMAT;
19 #endif
20 
21 static int my_max(int p, int q) { return (p > q) ? p : q; }
22 
23 static int my_min(int p, int q) { return (p < q) ? p : q; }
24 
25 
26 void BandFunctions(int l1, int u1, int l2, int u2)
27 {
28  int i, j;
29  BandMatrix BM1(20, l1, u1); BM1 = 999999.0;
30  for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
31  if (i - j <= l1 && i - j >= -u1) BM1(i, j) = 100 * i + j;
32 
33  BandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);
34 
35  BM2.ReSize(20, l2, u2); BM2 = 777777.0;
36  for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
37  if (i - j <= l2 && i - j >= -u2) BM2(i, j) = (100 * i + j) * 11;
38 
39  BandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
40  BMM = BM1 * BM2, BMN = -BM1;
41 
42  Matrix M1(20,20); M1 = 0;
43  for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
44  if (i - j <= l1 && i - j >= -u1) M1(i, j) = 100 * i + j;
45 
46  Matrix M2(20,20); M2 = 0;
47  for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
48  if (i - j <= l2 && i - j >= -u2) M2(i, j) = (100 * i + j) * 11;
49 
50  Matrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MM = M1 * M2, MN = -M1;
51  MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
52  Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);
53 
54  Matrix Test(7, 2);
55  Test(1,1) = BM1.BandWidth().Lower() - l1;
56  Test(1,2) = BM1.BandWidth().Upper() - u1;
57  Test(2,1) = BM2.BandWidth().Lower() - l2;
58  Test(2,2) = BM2.BandWidth().Upper() - u2;
59  Test(3,1) = BMA.BandWidth().Lower() - my_max(l1,l2);
60  Test(3,2) = BMA.BandWidth().Upper() - my_max(u1,u2);
61  Test(4,1) = BMS.BandWidth().Lower() - my_max(l1,l2);
62  Test(4,2) = BMS.BandWidth().Upper() - my_max(u1,u2);
63  Test(5,1) = BMSP.BandWidth().Lower() - my_min(l1,l2);
64  Test(5,2) = BMSP.BandWidth().Upper() - my_min(u1,u2);
65  Test(6,1) = BMM.BandWidth().Lower() - (l1 + l2);
66  Test(6,2) = BMM.BandWidth().Upper() - (u1 + u2);
67  Test(7,1) = BMN.BandWidth().Lower() - l1;
68  Test(7,2) = BMN.BandWidth().Upper() - u1;
69  Print(Test);
70 
71  // test element function
72  BandMatrix BM1E(20, l1, u1); BM1E = 999999.0;
73  for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
74  if (i - j <= l1 && i - j >= -u1) BM1E.element(i-1, j-1) = 100 * i + j;
75  BandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
76  for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
77  if (i - j <= l2 && i - j >= -u2)
78  BM2E.element(i-1, j-1) = (100 * i + j) * 11;
79  M1 = BM1E - BM1; Print(M1);
80  M2 = BM2E - BM2; Print(M2);
81 
82  // test element function with constant
83  BM1E = 444444.04; BM2E = 555555.0;
84  const BandMatrix BM1C = BM1, BM2C = BM2;
85  for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
86  if (i - j <= l1 && i - j >= -u1)
87  BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
88  for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
89  if (i - j <= l2 && i - j >= -u2)
90  BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
91  M1 = BM1E - BM1; Print(M1);
92  M2 = BM2E - BM2; Print(M2);
93 
94  // test subscript function with constant
95  BM1E = 444444.04; BM2E = 555555.0;
96  for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
97  if (i - j <= l1 && i - j >= -u1) BM1E(i, j) = BM1C(i, j);
98  for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
99  if (i - j <= l2 && i - j >= -u2) BM2E(i, j) = BM2C(i, j);
100  M1 = BM1E - BM1; Print(M1);
101  M2 = BM2E - BM2; Print(M2);
102 }
103 
104 void LowerBandFunctions(int l1, int l2)
105 {
106  int i, j;
107  LowerBandMatrix BM1(20, l1); BM1 = 999999.0;
108  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
109  if (i - j <= l1) BM1(i, j) = 100 * i + j;
110 
111  LowerBandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);
112 
113  BM2.ReSize(20, l2); BM2 = 777777.0;
114  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
115  if (i - j <= l2) BM2(i, j) = (100 * i + j) * 11;
116 
117  LowerBandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
118  BMM = BM1 * BM2, BMN = -BM1;
119 
120  Matrix M1(20,20); M1 = 0;
121  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
122  if (i - j <= l1) M1(i, j) = 100 * i + j;
123 
124  Matrix M2(20,20); M2 = 0;
125  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
126  if (i - j <= l2) M2(i, j) = (100 * i + j) * 11;
127 
128  Matrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MM = M1 * M2, MN = -M1;
129  MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
130  Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);
131 
132  Matrix Test(7, 2);
133  Test(1,1) = BM1.BandWidth().Lower() - l1;
134  Test(1,2) = BM1.BandWidth().Upper();
135  Test(2,1) = BM2.BandWidth().Lower() - l2;
136  Test(2,2) = BM2.BandWidth().Upper();
137  Test(3,1) = BMA.BandWidth().Lower() - my_max(l1,l2);
138  Test(3,2) = BMA.BandWidth().Upper();
139  Test(4,1) = BMS.BandWidth().Lower() - my_max(l1,l2);
140  Test(4,2) = BMS.BandWidth().Upper();
141  Test(5,1) = BMSP.BandWidth().Lower() - my_min(l1,l2);
142  Test(5,2) = BMSP.BandWidth().Upper();
143  Test(6,1) = BMM.BandWidth().Lower() - (l1 + l2);
144  Test(6,2) = BMM.BandWidth().Upper();
145  Test(7,1) = BMN.BandWidth().Lower() - l1;
146  Test(7,2) = BMN.BandWidth().Upper();
147  Print(Test);
148 
149  // test element function
150  LowerBandMatrix BM1E(20, l1); BM1E = 999999.0;
151  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
152  if (i - j <= l1) BM1E.element(i-1, j-1) = 100 * i + j;
153  LowerBandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
154  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
155  if (i - j <= l2) BM2E.element(i-1, j-1) = (100 * i + j) * 11;
156  M1 = BM1E - BM1; Print(M1);
157  M2 = BM2E - BM2; Print(M2);
158 
159  // test element function with constant
160  BM1E = 444444.04; BM2E = 555555.0;
161  const LowerBandMatrix BM1C = BM1, BM2C = BM2;
162  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
163  if (i - j <= l1) BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
164  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
165  if (i - j <= l2) BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
166  M1 = BM1E - BM1; Print(M1);
167  M2 = BM2E - BM2; Print(M2);
168 
169  // test subscript function with constant
170  BM1E = 444444.04; BM2E = 555555.0;
171  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
172  if (i - j <= l1) BM1E(i, j) = BM1C(i, j);
173  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
174  if (i - j <= l2) BM2E(i, j) = BM2C(i, j);
175  M1 = BM1E - BM1; Print(M1);
176  M2 = BM2E - BM2; Print(M2);
177 }
178 
179 void UpperBandFunctions(int u1, int u2)
180 {
181  int i, j;
182  UpperBandMatrix BM1(20, u1); BM1 = 999999.0;
183  for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
184  if (i - j >= -u1) BM1(i, j) = 100 * i + j;
185 
186  UpperBandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);
187 
188  BM2.ReSize(20, u2); BM2 = 777777.0;
189  for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
190  if (i - j >= -u2) BM2(i, j) = (100 * i + j) * 11;
191 
192  UpperBandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
193  BMM = BM1 * BM2, BMN = -BM1;
194 
195  Matrix M1(20,20); M1 = 0;
196  for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
197  if (i - j >= -u1) M1(i, j) = 100 * i + j;
198 
199  Matrix M2(20,20); M2 = 0;
200  for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
201  if (i - j >= -u2) M2(i, j) = (100 * i + j) * 11;
202 
203  Matrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MM = M1 * M2, MN = -M1;
204  MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
205  Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);
206 
207  Matrix Test(7, 2);
208  Test(1,1) = BM1.BandWidth().Lower();
209  Test(1,2) = BM1.BandWidth().Upper() - u1;
210  Test(2,1) = BM2.BandWidth().Lower();
211  Test(2,2) = BM2.BandWidth().Upper() - u2;
212  Test(3,1) = BMA.BandWidth().Lower();
213  Test(3,2) = BMA.BandWidth().Upper() - my_max(u1,u2);
214  Test(4,1) = BMS.BandWidth().Lower();
215  Test(4,2) = BMS.BandWidth().Upper() - my_max(u1,u2);
216  Test(5,1) = BMSP.BandWidth().Lower();
217  Test(5,2) = BMSP.BandWidth().Upper() - my_min(u1,u2);
218  Test(6,1) = BMM.BandWidth().Lower();
219  Test(6,2) = BMM.BandWidth().Upper() - (u1 + u2);
220  Test(7,1) = BMN.BandWidth().Lower();
221  Test(7,2) = BMN.BandWidth().Upper() - u1;
222  Print(Test);
223 
224  // test element function
225  UpperBandMatrix BM1E(20, u1); BM1E = 999999.0;
226  for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
227  if (i - j >= -u1) BM1E.element(i-1, j-1) = 100 * i + j;
228  UpperBandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
229  for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
230  if (i - j >= -u2) BM2E.element(i-1, j-1) = (100 * i + j) * 11;
231  M1 = BM1E - BM1; Print(M1);
232  M2 = BM2E - BM2; Print(M2);
233 
234  // test element function with constant
235  BM1E = 444444.04; BM2E = 555555.0;
236  const UpperBandMatrix BM1C = BM1, BM2C = BM2;
237  for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
238  if (i - j >= -u1) BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
239  for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
240  if (i - j >= -u2) BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
241  M1 = BM1E - BM1; Print(M1);
242  M2 = BM2E - BM2; Print(M2);
243 
244  // test subscript function with constant
245  BM1E = 444444.04; BM2E = 555555.0;
246  for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
247  if (i - j >= -u1) BM1E(i, j) = BM1C(i, j);
248  for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
249  if (i - j >= -u2) BM2E(i, j) = BM2C(i, j);
250  M1 = BM1E - BM1; Print(M1);
251  M2 = BM2E - BM2; Print(M2);
252 }
253 
254 void SymmetricBandFunctions(int l1, int l2)
255 {
256  int i, j;
257  SymmetricBandMatrix BM1(20, l1); BM1 = 999999.0;
258  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
259  if (i - j <= l1) BM1(i, j) = 100 * i + j;
260 
261  SymmetricBandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);
262 
263  BM2.ReSize(20, l2); BM2 = 777777.0;
264  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
265  if (i - j <= l2) BM2(i, j) = (100 * i + j) * 11;
266 
267  SymmetricBandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
268  BMN = -BM1;
269  BandMatrix BMM = BM1 * BM2;
270 
271  SymmetricMatrix M1(20); M1 = 0;
272  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
273  if (i - j <= l1) M1(i, j) = 100 * i + j;
274 
275  SymmetricMatrix M2(20); M2 = 0;
276  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
277  if (i - j <= l2) M2(i, j) = (100 * i + j) * 11;
278 
279  SymmetricMatrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MN = -M1;
280  Matrix MM = M1 * M2;
281  MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
282  Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);
283 
284  Matrix Test(7, 2);
285  Test(1,1) = BM1.BandWidth().Lower() - l1;
286  Test(1,2) = BM1.BandWidth().Upper() - l1;
287  Test(2,1) = BM2.BandWidth().Lower() - l2;
288  Test(2,2) = BM2.BandWidth().Upper() - l2;
289  Test(3,1) = BMA.BandWidth().Lower() - my_max(l1,l2);
290  Test(3,2) = BMA.BandWidth().Upper() - my_max(l1,l2);
291  Test(4,1) = BMS.BandWidth().Lower() - my_max(l1,l2);
292  Test(4,2) = BMS.BandWidth().Upper() - my_max(l1,l2);
293  Test(5,1) = BMSP.BandWidth().Lower() - my_min(l1,l2);
294  Test(5,2) = BMSP.BandWidth().Upper() - my_min(l1,l2);
295  Test(6,1) = BMM.BandWidth().Lower() - (l1 + l2);
296  Test(6,2) = BMM.BandWidth().Upper() - (l1 + l2);
297  Test(7,1) = BMN.BandWidth().Lower() - l1;
298  Test(7,2) = BMN.BandWidth().Upper() - l1;
299  Print(Test);
300 
301  // test element function
302  SymmetricBandMatrix BM1E(20, l1); BM1E = 999999.0;
303  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
304  if (i - j <= l1) BM1E.element(i-1, j-1) = 100 * i + j;
305  SymmetricBandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
306  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
307  if (i - j <= l2) BM2E.element(i-1, j-1) = (100 * i + j) * 11;
308  M1 = BM1E - BM1; Print(M1);
309  M2 = BM2E - BM2; Print(M2);
310 
311  // reverse subscripts
312  BM1E = 999999.0;
313  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
314  if (i - j <= l1) BM1E.element(j-1, i-1) = 100 * i + j;
315  BM2E = 777777.0;
316  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
317  if (i - j <= l2) BM2E.element(j-1, i-1) = (100 * i + j) * 11;
318  M1 = BM1E - BM1; Print(M1);
319  M2 = BM2E - BM2; Print(M2);
320 
321  // test element function with constant
322  BM1E = 444444.04; BM2E = 555555.0;
323  const SymmetricBandMatrix BM1C = BM1, BM2C = BM2;
324  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
325  if (i - j <= l1) BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
326  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
327  if (i - j <= l2) BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
328  M1 = BM1E - BM1; Print(M1);
329  M2 = BM2E - BM2; Print(M2);
330 
331  // reverse subscripts
332  BM1E = 444444.0; BM2E = 555555.0;
333  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
334  if (i - j <= l1) BM1E.element(j-1, i-1) = BM1C.element(j-1, i-1);
335  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
336  if (i - j <= l2) BM2E.element(j-1, i-1) = BM2C.element(j-1, i-1);
337  M1 = BM1E - BM1; Print(M1);
338  M2 = BM2E - BM2; Print(M2);
339 
340  // test subscript function with constant
341  BM1E = 444444.0; BM2E = 555555.0;
342  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
343  if (i - j <= l1) BM1E(i, j) = BM1C(i, j);
344  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
345  if (i - j <= l2) BM2E(i, j) = BM2C(i, j);
346  M1 = BM1E - BM1; Print(M1);
347  M2 = BM2E - BM2; Print(M2);
348 
349  // reverse subscripts
350  BM1E = 444444.0; BM2E = 555555.0;
351  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
352  if (i - j <= l1) BM1E(j, i) = BM1C(j, i);
353  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
354  if (i - j <= l2) BM2E(j, i) = BM2C(j, i);
355  M1 = BM1E - BM1; Print(M1);
356  M2 = BM2E - BM2; Print(M2);
357 
358  // partly reverse subscripts
359  BM1E = 444444.0; BM2E = 555555.0;
360  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
361  if (i - j <= l1) BM1E(j, i) = BM1C(i, j);
362  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
363  if (i - j <= l2) BM2E(j, i) = BM2C(i, j);
364  M1 = BM1E - BM1; Print(M1);
365  M2 = BM2E - BM2; Print(M2);
366 
367 }
368 
369 
370 
371 void trymath()
372 {
373 // cout << "\nSeventeenth test of Matrix package\n";
374 // cout << "\n";
375  Tracer et("Seventeenth test of Matrix package");
377 
378 
379  {
380  Tracer et1("Stage 1");
381  int i, j;
382  BandMatrix B(8,3,1);
383  for (i=1; i<=8; i++) for (j=-3; j<=1; j++)
384  { int k = i+j; if (k>0 && k<=8) B(i,k) = i + k/64.0; }
385 
386  IdentityMatrix I(8); BandMatrix B1; B1 = I;
387  UpperTriangularMatrix UT = I; Print(Matrix(B1-UT));
388  Print(Matrix(B * B - B * 2 + I - (B - I) * (B - I)));
389  Matrix A = B; BandMatrix C; C = B;
390  Print(Matrix(B * A - C * 2 + I - (A - I) * (B - I)));
391 
392  C.ReSize(8,2,2); C = 0.0; C.Inject(B);
393  Matrix X = A.t();
394  B1.ReSize(8,2,2); B1.Inject(X); Print(Matrix(C.t()-B1));
395 
396  Matrix BI = B.i(); A = A.i()-BI; Clean(A,0.000000001); Print(A);
397  BandLUMatrix BLU = B.t();
398  BI = BLU.i(); A = X.i()-BI; Clean(A,0.000000001); Print(A);
399  X.ReSize(1,1);
400  X(1,1) = BLU.LogDeterminant().Value()-B.LogDeterminant().Value();
401  Clean(X,0.000000001); Print(X);
402 
403  UpperBandMatrix U; U << B; LowerBandMatrix L; L << B;
404  DiagonalMatrix D; D << B;
405  Print( Matrix(L + (U - D - B)) );
406 
407 
408 
409  for (i=1; i<=8; i++) A.Column(i) << B.Column(i);
410  Print(Matrix(A-B));
411  }
412  {
413  Tracer et1("Stage 2");
414  BandMatrix A(7,2,2);
415  int i,j;
416  for (i=1; i<=7; i++) for (j=1; j<=7; j++)
417  {
418  int k=i-j; if (k<0) k = -k;
419  if (k==0) A(i,j)=6;
420  else if (k==1) A(i,j) = -4;
421  else if (k==2) A(i,j) = 1;
422  A(1,1) = A(7,7) = 5;
423  }
424  DiagonalMatrix D(7); D = 0.0; A = A - D;
425  BandLUMatrix B(A); Matrix M = A;
426  ColumnVector V(6);
427  V(1) = LogDeterminant(B).Value();
428  V(2) = LogDeterminant(A).Value();
429  V(3) = LogDeterminant(M).Value();
430  V(4) = Determinant(B);
431  V(5) = Determinant(A);
432  V(6) = Determinant(M);
433  V = V / 64 - 1; Clean(V,0.000000001); Print(V);
434  ColumnVector X(7);
435 
436  Real a[] = {1,2,3,4,5,6,7};
437  X << a;
438  M = (M.i()*X).t() - (B.i()*X).t() * 2.0 + (A.i()*X).t();
439  Clean(M,0.000000001); Print(M);
440 
441 
442  BandMatrix P(80,2,5); ColumnVector CX(80);
443  for (i=1; i<=80; i++) for (j=1; j<=80; j++)
444  { int d = i-j; if (d<=2 && d>=-5) P(i,j) = i + j/100.0; }
445  for (i=1; i<=80; i++) CX(i) = i*100.0;
446  Matrix MP = P;
447  ColumnVector V1 = P.i() * CX; ColumnVector V2 = MP.i() * CX;
448  V = V1 - V2; Clean(V,0.000000001); Print(V);
449 
450  V1 = P * V1; V2 = MP * V2; V = V1 - V2; Clean(V,0.000000001); Print(V);
451  RowVector XX(1);
452  XX = LogDeterminant(P).Value() / LogDeterminant(MP).Value() - 1.0;
453  Clean(XX,0.000000001); Print(XX);
454 
455  LowerBandMatrix LP(80,5);
456  for (i=1; i<=80; i++) for (j=1; j<=80; j++)
457  { int d = i-j; if (d<=5 && d>=0) LP(i,j) = i + j/100.0; }
458  MP = LP;
459  XX.ReSize(4);
460  XX(1) = LogDeterminant(LP).Value();
461  XX(2) = LogDeterminant(MP).Value();
462  V1 = LP.i() * CX; V2 = MP.i() * CX;
463  V = V1 - V2; Clean(V,0.000000001); Print(V);
464 
465  UpperBandMatrix UP(80,4);
466  for (i=1; i<=80; i++) for (j=1; j<=80; j++)
467  { int d = i-j; if (d<=0 && d>=-4) UP(i,j) = i + j/100.0; }
468  MP = UP;
469  XX(3) = LogDeterminant(UP).Value();
470  XX(4) = LogDeterminant(MP).Value();
471  V1 = UP.i() * CX; V2 = MP.i() * CX;
472  V = V1 - V2; Clean(V,0.000000001); Print(V);
473  XX = XX / SumAbsoluteValue(XX) - .25; Clean(XX,0.000000001); Print(XX);
474  }
475 
476  {
477  Tracer et1("Stage 3");
478  SymmetricBandMatrix SA(8,5);
479  int i,j;
480  for (i=1; i<=8; i++) for (j=1; j<=8; j++)
481  if (i-j<=5 && 0<=i-j) SA(i,j) =i + j/128.0;
482  DiagonalMatrix D(8); D = 10; SA = SA + D;
483 
484  Matrix MA1(8,8); Matrix MA2(8,8);
485  for (i=1; i<=8; i++)
486  { MA1.Column(i) << SA.Column(i); MA2.Row(i) << SA.Row(i); }
487  Print(Matrix(MA1-MA2));
488 
489  D = 10; SA = SA.t() + D; MA1 = MA1 + D;
490  Print(Matrix(MA1-SA));
491 
492  UpperBandMatrix UB(8,3); LowerBandMatrix LB(8,4);
493  D << SA; UB << SA; LB << SA;
494  SA = SA * 5.0; D = D * 5.0; LB = LB * 5.0; UB = UB * 5.0;
495  BandMatrix B = LB - D + UB - SA; Print(Matrix(B));
496 
497  SymmetricBandMatrix A(7,2); A = 100.0;
498  for (i=1; i<=7; i++) for (j=1; j<=7; j++)
499  {
500  int k=i-j;
501  if (k==0) A(i,j)=6;
502  else if (k==1) A(i,j) = -4;
503  else if (k==2) A(i,j) = 1;
504  A(1,1) = A(7,7) = 5;
505  }
506  BandLUMatrix C(A); Matrix M = A;
507  ColumnVector X(8);
508  X(1) = LogDeterminant(C).Value() - 64;
509  X(2) = LogDeterminant(A).Value() - 64;
510  X(3) = LogDeterminant(M).Value() - 64;
511  X(4) = SumSquare(M) - SumSquare(A);
512  X(5) = SumAbsoluteValue(M) - SumAbsoluteValue(A);
514  X(7) = Trace(M) - Trace(A);
515  X(8) = Sum(M) - Sum(A);
516  Clean(X,0.000000001); Print(X);
517 
518  Real a[] = {1,2,3,4,5,6,7};
519  X.ReSize(7);
520  X << a;
521  X = M.i()*X - C.i()*X * 2 + A.i()*X;
522  Clean(X,0.000000001); Print(X);
523 
524 
525  LB << A; UB << A; D << A;
526  BandMatrix XA = LB + (UB - D);
527  Print(Matrix(XA - A));
528 
529  for (i=1; i<=7; i++) for (j=1; j<=7; j++)
530  {
531  int k=i-j;
532  if (k==0) A(i,j)=6;
533  else if (k==1) A(i,j) = -4;
534  else if (k==2) A(i,j) = 1;
535  A(1,1) = A(7,7) = 5;
536  }
537  D = 1;
538 
539  M = LB.i() * LB - D; Clean(M,0.000000001); Print(M);
540  M = UB.i() * UB - D; Clean(M,0.000000001); Print(M);
541  M = XA.i() * XA - D; Clean(M,0.000000001); Print(M);
542  Matrix MUB = UB; Matrix MLB = LB;
543  M = LB.i() * UB - LB.i() * MUB; Clean(M,0.000000001); Print(M);
544  M = UB.i() * LB - UB.i() * MLB; Clean(M,0.000000001); Print(M);
545  M = LB.i() * UB - LB.i() * Matrix(UB); Clean(M,0.000000001); Print(M);
546  M = UB.i() * LB - UB.i() * Matrix(LB); Clean(M,0.000000001); Print(M);
547  }
548 
549  {
550  // some tests about adding and subtracting band matrices of different
551  // sizes - check bandwidth of results
552  Tracer et1("Stage 4");
553 
554  BandFunctions(9, 3, 9, 3); // equal
555  BandFunctions(4, 7, 4, 7); // equal
556  BandFunctions(9, 3, 5, 8); // neither < or >
557  BandFunctions(5, 8, 9, 3); // neither < or >
558  BandFunctions(9, 8, 5, 3); // >
559  BandFunctions(3, 5, 8, 9); // <
560 
561  LowerBandFunctions(9, 9); // equal
562  LowerBandFunctions(4, 4); // equal
563  LowerBandFunctions(9, 5); // >
564  LowerBandFunctions(3, 8); // <
565 
566  UpperBandFunctions(3, 3); // equal
567  UpperBandFunctions(7, 7); // equal
568  UpperBandFunctions(8, 3); // >
569  UpperBandFunctions(5, 9); // <
570 
571  SymmetricBandFunctions(9, 9); // equal
572  SymmetricBandFunctions(4, 4); // equal
573  SymmetricBandFunctions(9, 5); // >
574  SymmetricBandFunctions(3, 8); // <
575 
576  DiagonalMatrix D(6); D << 2 << 3 << 4.5 << 1.25 << 9.5 << -5;
577  BandMatrix BD = D;
578  UpperBandMatrix UBD; UBD = D;
579  LowerBandMatrix LBD; LBD = D;
580  SymmetricBandMatrix SBD = D;
581  Matrix X = BD - D; Print(X); X = UBD - D; Print(X);
582  X = LBD - D; Print(X); X = SBD - D; Print(X);
583  Matrix Test(9,2);
584  Test(1,1) = BD.BandWidth().Lower(); Test(1,2) = BD.BandWidth().Upper();
585  Test(2,1) = UBD.BandWidth().Lower(); Test(2,2) = UBD.BandWidth().Upper();
586  Test(3,1) = LBD.BandWidth().Lower(); Test(3,2) = LBD.BandWidth().Upper();
587  Test(4,1) = SBD.BandWidth().Lower(); Test(4,2) = SBD.BandWidth().Upper();
588 
589  IdentityMatrix I(10); I *= 5;
590  BD = I; UBD = I; LBD = I; SBD = I;
591  X = BD - I; Print(X); X = UBD - I; Print(X);
592  X = LBD - I; Print(X); X = SBD - I; Print(X);
593  Test(5,1) = BD.BandWidth().Lower(); Test(5,2) = BD.BandWidth().Upper();
594  Test(6,1) = UBD.BandWidth().Lower(); Test(6,2) = UBD.BandWidth().Upper();
595  Test(7,1) = LBD.BandWidth().Lower(); Test(7,2) = LBD.BandWidth().Upper();
596  Test(8,1) = SBD.BandWidth().Lower(); Test(8,2) = SBD.BandWidth().Upper();
597 
598  RowVector RV = D.AsRow(); I.ReSize(6); BandMatrix BI = I; I = 1;
599  BD = RV.AsDiagonal() + BI; X = BD - D - I; Print(X);
600  Test(9,1) = BD.BandWidth().Lower(); Test(9,2) = BD.BandWidth().Upper();
601 
602  Print(Test);
603  }
604 
605  {
606  // various element functions
607  Tracer et1("Stage 5");
608 
609  int i, j;
610  Matrix Count(1, 1); Count = 0; // for counting errors
611  Matrix M(20,30);
612  for (i = 1; i <= 20; ++i) for (j = 1; j <= 30; ++j)
613  M(i, j) = 100 * i + j;
614  const Matrix CM = M;
615  for (i = 1; i <= 20; ++i) for (j = 1; j <= 30; ++j)
616  {
617  if (M(i, j) != CM(i, j)) ++Count(1,1);
618  if (M(i, j) != CM.element(i-1, j-1)) ++Count(1,1);
619  if (M(i, j) != M.element(i-1, j-1)) ++Count(1,1);
620  }
621 
622  UpperTriangularMatrix U(20);
623  for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
624  U(i, j) = 100 * i + j;
625  const UpperTriangularMatrix CU = U;
626  for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
627  {
628  if (U(i, j) != CU(i, j)) ++Count(1,1);
629  if (U(i, j) != CU.element(i-1, j-1)) ++Count(1,1);
630  if (U(i, j) != U.element(i-1, j-1)) ++Count(1,1);
631  }
632 
633  LowerTriangularMatrix L(20);
634  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
635  L(i, j) = 100 * i + j;
636  const LowerTriangularMatrix CL = L;
637  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
638  {
639  if (L(i, j) != CL(i, j)) ++Count(1,1);
640  if (L(i, j) != CL.element(i-1, j-1)) ++Count(1,1);
641  if (L(i, j) != L.element(i-1, j-1)) ++Count(1,1);
642  }
643 
644  SymmetricMatrix S(20);
645  for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
646  S(i, j) = 100 * i + j;
647  const SymmetricMatrix CS = S;
648  for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
649  {
650  if (S(i, j) != CS(i, j)) ++Count(1,1);
651  if (S(i, j) != CS.element(i-1, j-1)) ++Count(1,1);
652  if (S(i, j) != S.element(i-1, j-1)) ++Count(1,1);
653  if (S(i, j) != S(j, i)) ++Count(1,1);
654  if (S(i, j) != CS(i, j)) ++Count(1,1);
655  if (S(i, j) != CS.element(i-1, j-1)) ++Count(1,1);
656  if (S(i, j) != S.element(i-1, j-1)) ++Count(1,1);
657  }
658 
659  DiagonalMatrix D(20);
660  for (i = 1; i <= 20; ++i) D(i) = 100 * i + i * i;
661  const DiagonalMatrix CD = D;
662  for (i = 1; i <= 20; ++i)
663  {
664  if (D(i, i) != CD(i, i)) ++Count(1,1);
665  if (D(i, i) != CD.element(i-1, i-1)) ++Count(1,1);
666  if (D(i, i) != D.element(i-1, i-1)) ++Count(1,1);
667  if (D(i, i) != D(i)) ++Count(1,1);
668  if (D(i) != CD(i)) ++Count(1,1);
669  if (D(i) != CD.element(i-1)) ++Count(1,1);
670  if (D(i) != D.element(i-1)) ++Count(1,1);
671  }
672 
673  RowVector R(20);
674  for (i = 1; i <= 20; ++i) R(i) = 100 * i + i * i;
675  const RowVector CR = R;
676  for (i = 1; i <= 20; ++i)
677  {
678  if (R(i) != CR(i)) ++Count(1,1);
679  if (R(i) != CR.element(i-1)) ++Count(1,1);
680  if (R(i) != R.element(i-1)) ++Count(1,1);
681  }
682 
683  ColumnVector C(20);
684  for (i = 1; i <= 20; ++i) C(i) = 100 * i + i * i;
685  const ColumnVector CC = C;
686  for (i = 1; i <= 20; ++i)
687  {
688  if (C(i) != CC(i)) ++Count(1,1);
689  if (C(i) != CC.element(i-1)) ++Count(1,1);
690  if (C(i) != C.element(i-1)) ++Count(1,1);
691  }
692 
693  Print(Count);
694 
695  }
696 
697  {
698  // resize to another matrix size
699  Tracer et1("Stage 6");
700 
701  Matrix A(20, 30); A = 3;
702  Matrix B(3, 4);
703  B.ReSize(A); B = 6; B -= 2 * A; Print(B);
704 
705  A.ReSize(25,25); A = 12;
706 
708  U.ReSize(A); U = 12; U << (U - A); Print(U);
709 
711  L.ReSize(U); L = 12; L << (L - A); Print(L);
712 
713  DiagonalMatrix D(5);
714  D.ReSize(U); D = 12; D << (D - A); Print(D);
715 
716  SymmetricMatrix S(5);
717  S.ReSize(U); S = 12; S << (S - A); Print(S);
718 
719  IdentityMatrix I(5);
720  I.ReSize(U); I = 12; D << (I - A); Print(D);
721 
722  A.ReSize(10, 1); A = 17;
723  ColumnVector C(5); C.ReSize(A); C = 17; C -= A; Print(C);
724 
725  A.ReSize(1, 10); A = 15;
726  RowVector R(5); R.ReSize(A); R = 15; R -= A; Print(R);
727 
728  }
729 
730  {
731  // generic matrix and identity matrix
732  Tracer et1("Stage 7");
733  IdentityMatrix I(5);
734  I *= 4;
735  GenericMatrix GM = I;
736  GM /= 2;
737  DiagonalMatrix D = GM;
738  Matrix A = GM + 10;
739  A -= 10;
740  A -= D;
741  Print(A);
742  }
743 
744  {
745  // SP and upper and lower triangular matrices
746  Tracer et1("Stage 8");
747  UpperTriangularMatrix UT(4);
748  UT << 3 << 7 << 3 << 9
749  << 5 << 2 << 6
750  << 8 << 0
751  << 4;
752  LowerTriangularMatrix LT; LT.ReSize(UT);
753  LT << 2
754  << 7 << 9
755  << 2 << 8 << 6
756  << 1 << 0 << 3 << 5;
757 
758  DiagonalMatrix D = SP(UT, LT);
759  DiagonalMatrix D1(4);
760  D1 << 6 << 45 << 48 << 20;
761  D -= D1; Print(D);
762  BandMatrix BM = SP(UT, LT);
763  Matrix X = BM - D1; Print(X);
764  RowVector RV(2);
765  RV(1) = BM.BandWidth().Lower();
766  RV(2) = BM.BandWidth().Upper();
767  Print(RV);
768  }
769 
770  {
771  // Helmert multiplies
772  Tracer et1("Stage 9");
773  MultWithCarry MCW;
774  int i, j;
775 
776  IdentityMatrix I(8);
777  Matrix X = I;
778  Matrix Y = Helmert_transpose(X);
779  Matrix H = Helmert(9); H -= Y.t(); Clean(H,0.000000001); Print(H);
780  Matrix Z = Helmert(Y) - I;
781  Clean(Z,0.000000001); Print(Z);
782 
783  Matrix A(9, 8);
784  for (i = 1; i <= 9; ++i) for (j = 1; j <= 8; ++j)
785  A(i, j) = Helmert_transpose(X.column(j), i);
786  A -= Y; Clean(A,0.000000001); Print(A);
787 
788  X = I;
789  Y = Helmert_transpose(X, true);
790  H = Helmert(8, true); H -= Y.t(); Clean(H,0.000000001); Print(H);
791  Z = Helmert(Y, true) - I;
792  Clean(Z,0.000000001); Print(Z);
793 
794  A.resize(8, 8);
795  for (i = 1; i <= 8; ++i) for (j = 1; j <= 8; ++j)
796  A(i, j) = Helmert_transpose(X.column(j), i, true);
797  A -= Y; Clean(A,0.000000001); Print(A);
798 
799 
800 
801  I.ReSize(9);
802  X = I;
803  Y = Helmert(X, true);
804  H = Helmert(9, true); H -= Y; Clean(H,0.000000001); Print(H);
805  Z = Helmert_transpose(Y, true) - I;
806  Clean(Z,0.000000001); Print(Z);
807 
808  A.ReSize(9, 9);
809  for (i = 1; i <= 9; ++i) A.Column(i) = Helmert(9, i, true);
810  A -= Y; Clean(A,0.000000001); Print(A);
811 
812  Y = Helmert(X);
813  A.ReSize(8, 9);
814  for (i = 1; i <= 9; ++i) A.Column(i) = Helmert(9, i);
815  A -= Y; Clean(A,0.000000001); Print(A);
816 
817  ColumnVector Twos(100); Twos = 2;
818  ColumnVector CV = Helmert(Twos); Clean(CV,0.000000001); Print(CV);
819 
820  X.resize(25,30);
821  FillWithValues(MCW, X);
822  Y = Helmert(X);
823  Z = Helmert(X,true).rows(1,24) - Y;
824  Clean(Z,0.000000001); Print(Z);
825  Z = Helmert(X,true).row(25) - X.sum_columns() / 5.0;
826  Clean(Z,0.000000001); Print(Z);
827 
828  I.resize(15);
829  X = I;
830  Z = Helmert_transpose(X, true) - Helmert(X, true).t();
831  Clean(Z,0.000000001); Print(Z);
832  I.resize(14); Y = I;
833  Z = Helmert(X) - Helmert_transpose(Y).t();
834  Clean(Z,0.000000001); Print(Z);
835 
836 
837 
838  }
839 
840 
841 
842 
843 
844 
845 // cout << "\nEnd of Seventeenth test\n";
846 }
847 
848 
LogAndSign LogDeterminant(const BaseMatrix &B)
Definition: newmat.h:2088
void ReSize(int m, int n, int b)
Definition: newmat.h:1224
Real & element(int, int)
Definition: newmat6.cpp:842
Real Sum(const BaseMatrix &B)
Definition: newmat.h:2107
ReturnMatrix Helmert_transpose(const ColumnVector &Y, bool full=false)
Definition: nm_misc.cpp:76
void ReSize(int m, int n, int b)
Definition: newmat.h:1185
LogAndSign LogDeterminant() const
Definition: newmat.h:342
void SymmetricBandFunctions(int l1, int l2)
Definition: tmth.cpp:254
void ReSize(int m)
Definition: newmat.h:786
void ReSize(int m)
Definition: newmat.h:1034
Real Value() const
Definition: newmat.h:64
void trymath()
Definition: tmth.cpp:371
Real & element(int, int)
Definition: newmat6.cpp:748
Real & element(int, int)
Definition: newmat6.cpp:732
ReturnMatrix Helmert(int n, bool full=false)
Definition: nm_misc.cpp:24
Real & element(int, int)
Definition: newmat6.cpp:716
GetSubMatrix rows(int, int) const
Definition: submat.cpp:58
void LowerBandFunctions(int l1, int l2)
Definition: tmth.cpp:104
SPMatrix SP(const BaseMatrix &, const BaseMatrix &)
Definition: newmat6.cpp:278
void resize(int n)
Definition: newmat4.cpp:335
virtual void ReSize(int m, int n)
Definition: newmat.h:662
void Inject(const GeneralMatrix &GM)
Definition: newmat.h:526
void ReSize(int m)
Definition: newmat.h:985
static int my_max(int p, int q)
Definition: tmth.cpp:21
double Real
Definition: include.h:307
Real & element(int, int)
Definition: newmat6.cpp:698
GetSubMatrix column(int) const
Definition: submat.cpp:68
Real & element(int, int)
Definition: newmat6.cpp:682
void ReSize(int m)
Definition: newmat.h:935
virtual void ReSize(int m, int n, int b)
Definition: newmat.h:1145
Real SumSquare(const BaseMatrix &B)
Definition: newmat.h:2095
void UpperBandFunctions(int u1, int u2)
Definition: tmth.cpp:179
virtual MatrixBandWidth BandWidth() const
Definition: newmat.h:389
Upper triangular band matrix.
Definition: newmat.h:1167
ReturnMatrix sum_columns() const
Definition: newmat8.cpp:805
Upper triangular matrix.
Definition: newmat.h:799
int Upper() const
Definition: newmat.h:216
Real & element(int, int)
Definition: newmat6.cpp:806
void Clean(Matrix &A, Real c)
Definition: tmt.cpp:128
Real & element(int, int)
Definition: newmat6.cpp:824
GetSubMatrix Column(int f) const
Definition: newmat.h:2152
FloatVector * d
Real & element(int, int)
Definition: newmat6.cpp:860
RowedMatrix AsRow() const
Definition: newmat.h:2141
Band matrix.
Definition: newmat.h:1096
TransposedMatrix t() const
Definition: newmat6.cpp:320
Real Trace(const BaseMatrix &B)
Definition: newmat.h:2100
Real MaximumAbsoluteValue(const BaseMatrix &B)
Definition: newmat.h:2111
void ReSize(int m)
Definition: newmat.h:881
void ReSize(int m, int b)
Definition: newmat.h:1289
static void PrintTrace()
Definition: myexcept.cpp:109
The usual rectangular matrix.
Definition: newmat.h:625
InvertedMatrix i() const
Definition: newmat6.cpp:329
static int my_min(int p, int q)
Definition: tmth.cpp:23
Real & element(int)
Definition: newmat6.cpp:778
Real SumAbsoluteValue(const BaseMatrix &B)
Definition: newmat.h:2103
FloatVector FloatVector * a
LU decomposition of a band matrix.
Definition: newmat.h:1306
Diagonal matrix.
Definition: newmat.h:896
Lower triangular matrix.
Definition: newmat.h:848
GetSubMatrix Row(int f) const
Definition: newmat.h:2150
Symmetric band matrix.
Definition: newmat.h:1245
Lower triangular band matrix.
Definition: newmat.h:1206
virtual void resize(int, int)
Definition: newmat4.cpp:289
DiagedMatrix AsDiagonal() const
Definition: newmat.h:2143
Real Determinant(const BaseMatrix &B)
Definition: newmat.h:2092
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
void BandFunctions(int l1, int u1, int l2, int u2)
Definition: tmth.cpp:26
void ReSize(int m)
Definition: newmat.h:833
GetSubMatrix row(int) const
Definition: submat.cpp:49
A matrix which can be of any GeneralMatrix type.
Definition: newmat.h:1397
Real & element(int)
Definition: newmat6.cpp:792
Identity matrix.
Definition: newmat.h:1350
void FillWithValues(MultWithCarry &MWC, Matrix &M)
Definition: tmt.cpp:190
void ReSize(int n)
Definition: newmat.h:1381
Symmetric matrix.
Definition: newmat.h:753


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