Planet
navi homePPSaboutscreenshotsdownloaddevelopmentforum

source: orxonox.OLD/orxonox/trunk/src/lib/newmat/tmth.cpp @ 4569

Last change on this file since 4569 was 4565, checked in by patrick, 19 years ago

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