1 /// \ingroup newmat
2 ///@{
3
4 /// \file tmth.cpp
5 /// Part of matrix library test program.
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
my_max(int p,int q)21 static int my_max(int p, int q) { return (p > q) ? p : q; }
22
my_min(int p,int q)23 static int my_min(int p, int q) { return (p < q) ? p : q; }
24
25
BandFunctions(int l1,int u1,int l2,int u2)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
LowerBandFunctions(int l1,int l2)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
UpperBandFunctions(int u1,int u2)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
SymmetricBandFunctions(int l1,int l2)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
trymath()371 void trymath()
372 {
373 // cout << "\nSeventeenth test of Matrix package\n";
374 // cout << "\n";
375 Tracer et("Seventeenth test of Matrix package");
376 Tracer::PrintTrace();
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);
513 X(6) = MaximumAbsoluteValue(M) - MaximumAbsoluteValue(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
707 UpperTriangularMatrix U(5);
708 U.ReSize(A); U = 12; U << (U - A); Print(U);
709
710 LowerTriangularMatrix L(5);
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
849 ///*}
850