1 | /// \ingroup newmat |
---|
2 | ///@{ |
---|
3 | |
---|
4 | /// \file newmat8.cpp |
---|
5 | /// LU transform, scalar functions of matrices. |
---|
6 | |
---|
7 | // Copyright (C) 1991,2,3,4,8: R B Davies |
---|
8 | |
---|
9 | #define WANT_MATH |
---|
10 | |
---|
11 | #include "include.h" |
---|
12 | |
---|
13 | #include "newmat.h" |
---|
14 | #include "newmatrc.h" |
---|
15 | #include "precisio.h" |
---|
16 | |
---|
17 | #ifdef use_namespace |
---|
18 | namespace NEWMAT { |
---|
19 | #endif |
---|
20 | |
---|
21 | |
---|
22 | #ifdef DO_REPORT |
---|
23 | #define REPORT { static ExeCounter ExeCount(__LINE__,8); ++ExeCount; } |
---|
24 | #else |
---|
25 | #define REPORT {} |
---|
26 | #endif |
---|
27 | |
---|
28 | |
---|
29 | /************************** LU transformation ****************************/ |
---|
30 | |
---|
31 | void CroutMatrix::ludcmp() |
---|
32 | // LU decomposition from Golub & Van Loan, algorithm 3.4.1, (the "outer |
---|
33 | // product" version). |
---|
34 | // This replaces the code derived from Numerical Recipes in C in previous |
---|
35 | // versions of newmat and being row oriented runs much faster with large |
---|
36 | // matrices. |
---|
37 | { |
---|
38 | REPORT |
---|
39 | Tracer tr( "Crout(ludcmp)" ); sing = false; |
---|
40 | Real* akk = store; // runs down diagonal |
---|
41 | |
---|
42 | Real big = fabs(*akk); int mu = 0; Real* ai = akk; int k; |
---|
43 | |
---|
44 | for (k = 1; k < nrows_val; k++) |
---|
45 | { |
---|
46 | ai += nrows_val; const Real trybig = fabs(*ai); |
---|
47 | if (big < trybig) { big = trybig; mu = k; } |
---|
48 | } |
---|
49 | |
---|
50 | |
---|
51 | if (nrows_val) for (k = 0;;) |
---|
52 | { |
---|
53 | /* |
---|
54 | int mu1; |
---|
55 | { |
---|
56 | Real big = fabs(*akk); mu1 = k; Real* ai = akk; int i; |
---|
57 | |
---|
58 | for (i = k+1; i < nrows_val; i++) |
---|
59 | { |
---|
60 | ai += nrows_val; const Real trybig = fabs(*ai); |
---|
61 | if (big < trybig) { big = trybig; mu1 = i; } |
---|
62 | } |
---|
63 | } |
---|
64 | if (mu1 != mu) cout << k << " " << mu << " " << mu1 << endl; |
---|
65 | */ |
---|
66 | |
---|
67 | indx[k] = mu; |
---|
68 | |
---|
69 | if (mu != k) //row swap |
---|
70 | { |
---|
71 | Real* a1 = store + nrows_val * k; |
---|
72 | Real* a2 = store + nrows_val * mu; d = !d; |
---|
73 | int j = nrows_val; |
---|
74 | while (j--) { const Real temp = *a1; *a1++ = *a2; *a2++ = temp; } |
---|
75 | } |
---|
76 | |
---|
77 | Real diag = *akk; big = 0; mu = k + 1; |
---|
78 | if (diag != 0) |
---|
79 | { |
---|
80 | ai = akk; int i = nrows_val - k - 1; |
---|
81 | while (i--) |
---|
82 | { |
---|
83 | ai += nrows_val; Real* al = ai; |
---|
84 | Real mult = *al / diag; *al = mult; |
---|
85 | int l = nrows_val - k - 1; Real* aj = akk; |
---|
86 | // work out the next pivot as part of this loop |
---|
87 | // this saves a column operation |
---|
88 | if (l-- != 0) |
---|
89 | { |
---|
90 | *(++al) -= (mult * *(++aj)); |
---|
91 | const Real trybig = fabs(*al); |
---|
92 | if (big < trybig) { big = trybig; mu = nrows_val - i - 1; } |
---|
93 | while (l--) *(++al) -= (mult * *(++aj)); |
---|
94 | } |
---|
95 | } |
---|
96 | } |
---|
97 | else sing = true; |
---|
98 | if (++k == nrows_val) break; // so next line won't overflow |
---|
99 | akk += nrows_val + 1; |
---|
100 | } |
---|
101 | } |
---|
102 | |
---|
103 | void CroutMatrix::lubksb(Real* B, int mini) |
---|
104 | { |
---|
105 | REPORT |
---|
106 | // this has been adapted from Numerical Recipes in C. The code has been |
---|
107 | // substantially streamlined, so I do not think much of the original |
---|
108 | // copyright remains. However there is not much opportunity for |
---|
109 | // variation in the code, so it is still similar to the NR code. |
---|
110 | // I follow the NR code in skipping over initial zeros in the B vector. |
---|
111 | |
---|
112 | Tracer tr("Crout(lubksb)"); |
---|
113 | if (sing) Throw(SingularException(*this)); |
---|
114 | int i, j, ii = nrows_val; // ii initialised : B might be all zeros |
---|
115 | |
---|
116 | |
---|
117 | // scan for first non-zero in B |
---|
118 | for (i = 0; i < nrows_val; i++) |
---|
119 | { |
---|
120 | int ip = indx[i]; Real temp = B[ip]; B[ip] = B[i]; B[i] = temp; |
---|
121 | if (temp != 0.0) { ii = i; break; } |
---|
122 | } |
---|
123 | |
---|
124 | Real* bi; Real* ai; |
---|
125 | i = ii + 1; |
---|
126 | |
---|
127 | if (i < nrows_val) |
---|
128 | { |
---|
129 | bi = B + ii; ai = store + ii + i * nrows_val; |
---|
130 | for (;;) |
---|
131 | { |
---|
132 | int ip = indx[i]; Real sum = B[ip]; B[ip] = B[i]; |
---|
133 | Real* aij = ai; Real* bj = bi; j = i - ii; |
---|
134 | while (j--) sum -= *aij++ * *bj++; |
---|
135 | B[i] = sum; |
---|
136 | if (++i == nrows_val) break; |
---|
137 | ai += nrows_val; |
---|
138 | } |
---|
139 | } |
---|
140 | |
---|
141 | ai = store + nrows_val * nrows_val; |
---|
142 | |
---|
143 | for (i = nrows_val - 1; i >= mini; i--) |
---|
144 | { |
---|
145 | Real* bj = B+i; ai -= nrows_val; Real* ajx = ai+i; |
---|
146 | Real sum = *bj; Real diag = *ajx; |
---|
147 | j = nrows_val - i; while(--j) sum -= *(++ajx) * *(++bj); |
---|
148 | B[i] = sum / diag; |
---|
149 | } |
---|
150 | } |
---|
151 | |
---|
152 | /****************************** scalar functions ****************************/ |
---|
153 | |
---|
154 | inline Real square(Real x) { return x*x; } |
---|
155 | |
---|
156 | Real GeneralMatrix::sum_square() const |
---|
157 | { |
---|
158 | REPORT |
---|
159 | Real sum = 0.0; int i = storage; Real* s = store; |
---|
160 | while (i--) sum += square(*s++); |
---|
161 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
162 | } |
---|
163 | |
---|
164 | Real GeneralMatrix::sum_absolute_value() const |
---|
165 | { |
---|
166 | REPORT |
---|
167 | Real sum = 0.0; int i = storage; Real* s = store; |
---|
168 | while (i--) sum += fabs(*s++); |
---|
169 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
170 | } |
---|
171 | |
---|
172 | Real GeneralMatrix::sum() const |
---|
173 | { |
---|
174 | REPORT |
---|
175 | Real sm = 0.0; int i = storage; Real* s = store; |
---|
176 | while (i--) sm += *s++; |
---|
177 | ((GeneralMatrix&)*this).tDelete(); return sm; |
---|
178 | } |
---|
179 | |
---|
180 | // maxima and minima |
---|
181 | |
---|
182 | // There are three sets of routines |
---|
183 | // maximum_absolute_value, minimum_absolute_value, maximum, minimum |
---|
184 | // ... these find just the maxima and minima |
---|
185 | // maximum_absolute_value1, minimum_absolute_value1, maximum1, minimum1 |
---|
186 | // ... these find the maxima and minima and their locations in a |
---|
187 | // one dimensional object |
---|
188 | // maximum_absolute_value2, minimum_absolute_value2, maximum2, minimum2 |
---|
189 | // ... these find the maxima and minima and their locations in a |
---|
190 | // two dimensional object |
---|
191 | |
---|
192 | // If the matrix has no values throw an exception |
---|
193 | |
---|
194 | // If we do not want the location find the maximum or minimum on the |
---|
195 | // array stored by GeneralMatrix |
---|
196 | // This won't work for BandMatrices. We call ClearCorner for |
---|
197 | // maximum_absolute_value but for the others use the absolute_minimum_value2 |
---|
198 | // version and discard the location. |
---|
199 | |
---|
200 | // For one dimensional objects, when we want the location of the |
---|
201 | // maximum or minimum, work with the array stored by GeneralMatrix |
---|
202 | |
---|
203 | // For two dimensional objects where we want the location of the maximum or |
---|
204 | // minimum proceed as follows: |
---|
205 | |
---|
206 | // For rectangular matrices use the array stored by GeneralMatrix and |
---|
207 | // deduce the location from the location in the GeneralMatrix |
---|
208 | |
---|
209 | // For other two dimensional matrices use the Matrix Row routine to find the |
---|
210 | // maximum or minimum for each row. |
---|
211 | |
---|
212 | static void NullMatrixError(const GeneralMatrix* gm) |
---|
213 | { |
---|
214 | ((GeneralMatrix&)*gm).tDelete(); |
---|
215 | Throw(ProgramException("Maximum or minimum of null matrix")); |
---|
216 | } |
---|
217 | |
---|
218 | Real GeneralMatrix::maximum_absolute_value() const |
---|
219 | { |
---|
220 | REPORT |
---|
221 | if (storage == 0) NullMatrixError(this); |
---|
222 | Real maxval = 0.0; int l = storage; Real* s = store; |
---|
223 | while (l--) { Real a = fabs(*s++); if (maxval < a) maxval = a; } |
---|
224 | ((GeneralMatrix&)*this).tDelete(); return maxval; |
---|
225 | } |
---|
226 | |
---|
227 | Real GeneralMatrix::maximum_absolute_value1(int& i) const |
---|
228 | { |
---|
229 | REPORT |
---|
230 | if (storage == 0) NullMatrixError(this); |
---|
231 | Real maxval = 0.0; int l = storage; Real* s = store; int li = storage; |
---|
232 | while (l--) |
---|
233 | { Real a = fabs(*s++); if (maxval <= a) { maxval = a; li = l; } } |
---|
234 | i = storage - li; |
---|
235 | ((GeneralMatrix&)*this).tDelete(); return maxval; |
---|
236 | } |
---|
237 | |
---|
238 | Real GeneralMatrix::minimum_absolute_value() const |
---|
239 | { |
---|
240 | REPORT |
---|
241 | if (storage == 0) NullMatrixError(this); |
---|
242 | int l = storage - 1; Real* s = store; Real minval = fabs(*s++); |
---|
243 | while (l--) { Real a = fabs(*s++); if (minval > a) minval = a; } |
---|
244 | ((GeneralMatrix&)*this).tDelete(); return minval; |
---|
245 | } |
---|
246 | |
---|
247 | Real GeneralMatrix::minimum_absolute_value1(int& i) const |
---|
248 | { |
---|
249 | REPORT |
---|
250 | if (storage == 0) NullMatrixError(this); |
---|
251 | int l = storage - 1; Real* s = store; Real minval = fabs(*s++); int li = l; |
---|
252 | while (l--) |
---|
253 | { Real a = fabs(*s++); if (minval >= a) { minval = a; li = l; } } |
---|
254 | i = storage - li; |
---|
255 | ((GeneralMatrix&)*this).tDelete(); return minval; |
---|
256 | } |
---|
257 | |
---|
258 | Real GeneralMatrix::maximum() const |
---|
259 | { |
---|
260 | REPORT |
---|
261 | if (storage == 0) NullMatrixError(this); |
---|
262 | int l = storage - 1; Real* s = store; Real maxval = *s++; |
---|
263 | while (l--) { Real a = *s++; if (maxval < a) maxval = a; } |
---|
264 | ((GeneralMatrix&)*this).tDelete(); return maxval; |
---|
265 | } |
---|
266 | |
---|
267 | Real GeneralMatrix::maximum1(int& i) const |
---|
268 | { |
---|
269 | REPORT |
---|
270 | if (storage == 0) NullMatrixError(this); |
---|
271 | int l = storage - 1; Real* s = store; Real maxval = *s++; int li = l; |
---|
272 | while (l--) { Real a = *s++; if (maxval <= a) { maxval = a; li = l; } } |
---|
273 | i = storage - li; |
---|
274 | ((GeneralMatrix&)*this).tDelete(); return maxval; |
---|
275 | } |
---|
276 | |
---|
277 | Real GeneralMatrix::minimum() const |
---|
278 | { |
---|
279 | REPORT |
---|
280 | if (storage == 0) NullMatrixError(this); |
---|
281 | int l = storage - 1; Real* s = store; Real minval = *s++; |
---|
282 | while (l--) { Real a = *s++; if (minval > a) minval = a; } |
---|
283 | ((GeneralMatrix&)*this).tDelete(); return minval; |
---|
284 | } |
---|
285 | |
---|
286 | Real GeneralMatrix::minimum1(int& i) const |
---|
287 | { |
---|
288 | REPORT |
---|
289 | if (storage == 0) NullMatrixError(this); |
---|
290 | int l = storage - 1; Real* s = store; Real minval = *s++; int li = l; |
---|
291 | while (l--) { Real a = *s++; if (minval >= a) { minval = a; li = l; } } |
---|
292 | i = storage - li; |
---|
293 | ((GeneralMatrix&)*this).tDelete(); return minval; |
---|
294 | } |
---|
295 | |
---|
296 | Real GeneralMatrix::maximum_absolute_value2(int& i, int& j) const |
---|
297 | { |
---|
298 | REPORT |
---|
299 | if (storage == 0) NullMatrixError(this); |
---|
300 | Real maxval = 0.0; int nr = Nrows(); |
---|
301 | MatrixRow mr((GeneralMatrix*)this, LoadOnEntry+DirectPart); |
---|
302 | for (int r = 1; r <= nr; r++) |
---|
303 | { |
---|
304 | int c; maxval = mr.MaximumAbsoluteValue1(maxval, c); |
---|
305 | if (c > 0) { i = r; j = c; } |
---|
306 | mr.Next(); |
---|
307 | } |
---|
308 | ((GeneralMatrix&)*this).tDelete(); return maxval; |
---|
309 | } |
---|
310 | |
---|
311 | Real GeneralMatrix::minimum_absolute_value2(int& i, int& j) const |
---|
312 | { |
---|
313 | REPORT |
---|
314 | if (storage == 0) NullMatrixError(this); |
---|
315 | Real minval = FloatingPointPrecision::Maximum(); int nr = Nrows(); |
---|
316 | MatrixRow mr((GeneralMatrix*)this, LoadOnEntry+DirectPart); |
---|
317 | for (int r = 1; r <= nr; r++) |
---|
318 | { |
---|
319 | int c; minval = mr.MinimumAbsoluteValue1(minval, c); |
---|
320 | if (c > 0) { i = r; j = c; } |
---|
321 | mr.Next(); |
---|
322 | } |
---|
323 | ((GeneralMatrix&)*this).tDelete(); return minval; |
---|
324 | } |
---|
325 | |
---|
326 | Real GeneralMatrix::maximum2(int& i, int& j) const |
---|
327 | { |
---|
328 | REPORT |
---|
329 | if (storage == 0) NullMatrixError(this); |
---|
330 | Real maxval = -FloatingPointPrecision::Maximum(); int nr = Nrows(); |
---|
331 | MatrixRow mr((GeneralMatrix*)this, LoadOnEntry+DirectPart); |
---|
332 | for (int r = 1; r <= nr; r++) |
---|
333 | { |
---|
334 | int c; maxval = mr.Maximum1(maxval, c); |
---|
335 | if (c > 0) { i = r; j = c; } |
---|
336 | mr.Next(); |
---|
337 | } |
---|
338 | ((GeneralMatrix&)*this).tDelete(); return maxval; |
---|
339 | } |
---|
340 | |
---|
341 | Real GeneralMatrix::minimum2(int& i, int& j) const |
---|
342 | { |
---|
343 | REPORT |
---|
344 | if (storage == 0) NullMatrixError(this); |
---|
345 | Real minval = FloatingPointPrecision::Maximum(); int nr = Nrows(); |
---|
346 | MatrixRow mr((GeneralMatrix*)this, LoadOnEntry+DirectPart); |
---|
347 | for (int r = 1; r <= nr; r++) |
---|
348 | { |
---|
349 | int c; minval = mr.Minimum1(minval, c); |
---|
350 | if (c > 0) { i = r; j = c; } |
---|
351 | mr.Next(); |
---|
352 | } |
---|
353 | ((GeneralMatrix&)*this).tDelete(); return minval; |
---|
354 | } |
---|
355 | |
---|
356 | Real Matrix::maximum_absolute_value2(int& i, int& j) const |
---|
357 | { |
---|
358 | REPORT |
---|
359 | int k; Real m = GeneralMatrix::maximum_absolute_value1(k); k--; |
---|
360 | i = k / Ncols(); j = k - i * Ncols(); i++; j++; |
---|
361 | return m; |
---|
362 | } |
---|
363 | |
---|
364 | Real Matrix::minimum_absolute_value2(int& i, int& j) const |
---|
365 | { |
---|
366 | REPORT |
---|
367 | int k; Real m = GeneralMatrix::minimum_absolute_value1(k); k--; |
---|
368 | i = k / Ncols(); j = k - i * Ncols(); i++; j++; |
---|
369 | return m; |
---|
370 | } |
---|
371 | |
---|
372 | Real Matrix::maximum2(int& i, int& j) const |
---|
373 | { |
---|
374 | REPORT |
---|
375 | int k; Real m = GeneralMatrix::maximum1(k); k--; |
---|
376 | i = k / Ncols(); j = k - i * Ncols(); i++; j++; |
---|
377 | return m; |
---|
378 | } |
---|
379 | |
---|
380 | Real Matrix::minimum2(int& i, int& j) const |
---|
381 | { |
---|
382 | REPORT |
---|
383 | int k; Real m = GeneralMatrix::minimum1(k); k--; |
---|
384 | i = k / Ncols(); j = k - i * Ncols(); i++; j++; |
---|
385 | return m; |
---|
386 | } |
---|
387 | |
---|
388 | Real SymmetricMatrix::sum_square() const |
---|
389 | { |
---|
390 | REPORT |
---|
391 | Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows_val; |
---|
392 | for (int i = 0; i<nr; i++) |
---|
393 | { |
---|
394 | int j = i; |
---|
395 | while (j--) sum2 += square(*s++); |
---|
396 | sum1 += square(*s++); |
---|
397 | } |
---|
398 | ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2; |
---|
399 | } |
---|
400 | |
---|
401 | Real SymmetricMatrix::sum_absolute_value() const |
---|
402 | { |
---|
403 | REPORT |
---|
404 | Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows_val; |
---|
405 | for (int i = 0; i<nr; i++) |
---|
406 | { |
---|
407 | int j = i; |
---|
408 | while (j--) sum2 += fabs(*s++); |
---|
409 | sum1 += fabs(*s++); |
---|
410 | } |
---|
411 | ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2; |
---|
412 | } |
---|
413 | |
---|
414 | Real IdentityMatrix::sum_absolute_value() const |
---|
415 | { REPORT return fabs(trace()); } // no need to do tDelete? |
---|
416 | |
---|
417 | Real SymmetricMatrix::sum() const |
---|
418 | { |
---|
419 | REPORT |
---|
420 | Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows_val; |
---|
421 | for (int i = 0; i<nr; i++) |
---|
422 | { |
---|
423 | int j = i; |
---|
424 | while (j--) sum2 += *s++; |
---|
425 | sum1 += *s++; |
---|
426 | } |
---|
427 | ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2; |
---|
428 | } |
---|
429 | |
---|
430 | Real IdentityMatrix::sum_square() const |
---|
431 | { |
---|
432 | Real sum = *store * *store * nrows_val; |
---|
433 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
434 | } |
---|
435 | |
---|
436 | |
---|
437 | Real BaseMatrix::sum_square() const |
---|
438 | { |
---|
439 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
440 | Real s = gm->sum_square(); return s; |
---|
441 | } |
---|
442 | |
---|
443 | Real BaseMatrix::norm_Frobenius() const |
---|
444 | { REPORT return sqrt(sum_square()); } |
---|
445 | |
---|
446 | Real BaseMatrix::sum_absolute_value() const |
---|
447 | { |
---|
448 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
449 | Real s = gm->sum_absolute_value(); return s; |
---|
450 | } |
---|
451 | |
---|
452 | Real BaseMatrix::sum() const |
---|
453 | { |
---|
454 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
455 | Real s = gm->sum(); return s; |
---|
456 | } |
---|
457 | |
---|
458 | Real BaseMatrix::maximum_absolute_value() const |
---|
459 | { |
---|
460 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
461 | Real s = gm->maximum_absolute_value(); return s; |
---|
462 | } |
---|
463 | |
---|
464 | Real BaseMatrix::maximum_absolute_value1(int& i) const |
---|
465 | { |
---|
466 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
467 | Real s = gm->maximum_absolute_value1(i); return s; |
---|
468 | } |
---|
469 | |
---|
470 | Real BaseMatrix::maximum_absolute_value2(int& i, int& j) const |
---|
471 | { |
---|
472 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
473 | Real s = gm->maximum_absolute_value2(i, j); return s; |
---|
474 | } |
---|
475 | |
---|
476 | Real BaseMatrix::minimum_absolute_value() const |
---|
477 | { |
---|
478 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
479 | Real s = gm->minimum_absolute_value(); return s; |
---|
480 | } |
---|
481 | |
---|
482 | Real BaseMatrix::minimum_absolute_value1(int& i) const |
---|
483 | { |
---|
484 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
485 | Real s = gm->minimum_absolute_value1(i); return s; |
---|
486 | } |
---|
487 | |
---|
488 | Real BaseMatrix::minimum_absolute_value2(int& i, int& j) const |
---|
489 | { |
---|
490 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
491 | Real s = gm->minimum_absolute_value2(i, j); return s; |
---|
492 | } |
---|
493 | |
---|
494 | Real BaseMatrix::maximum() const |
---|
495 | { |
---|
496 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
497 | Real s = gm->maximum(); return s; |
---|
498 | } |
---|
499 | |
---|
500 | Real BaseMatrix::maximum1(int& i) const |
---|
501 | { |
---|
502 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
503 | Real s = gm->maximum1(i); return s; |
---|
504 | } |
---|
505 | |
---|
506 | Real BaseMatrix::maximum2(int& i, int& j) const |
---|
507 | { |
---|
508 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
509 | Real s = gm->maximum2(i, j); return s; |
---|
510 | } |
---|
511 | |
---|
512 | Real BaseMatrix::minimum() const |
---|
513 | { |
---|
514 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
515 | Real s = gm->minimum(); return s; |
---|
516 | } |
---|
517 | |
---|
518 | Real BaseMatrix::minimum1(int& i) const |
---|
519 | { |
---|
520 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
521 | Real s = gm->minimum1(i); return s; |
---|
522 | } |
---|
523 | |
---|
524 | Real BaseMatrix::minimum2(int& i, int& j) const |
---|
525 | { |
---|
526 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
527 | Real s = gm->minimum2(i, j); return s; |
---|
528 | } |
---|
529 | |
---|
530 | Real dotproduct(const Matrix& A, const Matrix& B) |
---|
531 | { |
---|
532 | REPORT |
---|
533 | int n = A.storage; |
---|
534 | if (n != B.storage) |
---|
535 | { |
---|
536 | Tracer tr("dotproduct"); |
---|
537 | Throw(IncompatibleDimensionsException(A,B)); |
---|
538 | } |
---|
539 | Real sum = 0.0; Real* a = A.store; Real* b = B.store; |
---|
540 | while (n--) sum += *a++ * *b++; |
---|
541 | return sum; |
---|
542 | } |
---|
543 | |
---|
544 | Real Matrix::trace() const |
---|
545 | { |
---|
546 | REPORT |
---|
547 | Tracer tr("trace"); |
---|
548 | int i = nrows_val; int d = i+1; |
---|
549 | if (i != ncols_val) Throw(NotSquareException(*this)); |
---|
550 | Real sum = 0.0; Real* s = store; |
---|
551 | // while (i--) { sum += *s; s += d; } |
---|
552 | if (i) for (;;) { sum += *s; if (!(--i)) break; s += d; } |
---|
553 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
554 | } |
---|
555 | |
---|
556 | Real DiagonalMatrix::trace() const |
---|
557 | { |
---|
558 | REPORT |
---|
559 | int i = nrows_val; Real sum = 0.0; Real* s = store; |
---|
560 | while (i--) sum += *s++; |
---|
561 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
562 | } |
---|
563 | |
---|
564 | Real SymmetricMatrix::trace() const |
---|
565 | { |
---|
566 | REPORT |
---|
567 | int i = nrows_val; Real sum = 0.0; Real* s = store; int j = 2; |
---|
568 | // while (i--) { sum += *s; s += j++; } |
---|
569 | if (i) for (;;) { sum += *s; if (!(--i)) break; s += j++; } |
---|
570 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
571 | } |
---|
572 | |
---|
573 | Real LowerTriangularMatrix::trace() const |
---|
574 | { |
---|
575 | REPORT |
---|
576 | int i = nrows_val; Real sum = 0.0; Real* s = store; int j = 2; |
---|
577 | // while (i--) { sum += *s; s += j++; } |
---|
578 | if (i) for (;;) { sum += *s; if (!(--i)) break; s += j++; } |
---|
579 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
580 | } |
---|
581 | |
---|
582 | Real UpperTriangularMatrix::trace() const |
---|
583 | { |
---|
584 | REPORT |
---|
585 | int i = nrows_val; Real sum = 0.0; Real* s = store; |
---|
586 | while (i) { sum += *s; s += i--; } // won t cause a problem |
---|
587 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
588 | } |
---|
589 | |
---|
590 | Real BandMatrix::trace() const |
---|
591 | { |
---|
592 | REPORT |
---|
593 | int i = nrows_val; int w = lower_val+upper_val+1; |
---|
594 | Real sum = 0.0; Real* s = store+lower_val; |
---|
595 | // while (i--) { sum += *s; s += w; } |
---|
596 | if (i) for (;;) { sum += *s; if (!(--i)) break; s += w; } |
---|
597 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
598 | } |
---|
599 | |
---|
600 | Real SymmetricBandMatrix::trace() const |
---|
601 | { |
---|
602 | REPORT |
---|
603 | int i = nrows_val; int w = lower_val+1; |
---|
604 | Real sum = 0.0; Real* s = store+lower_val; |
---|
605 | // while (i--) { sum += *s; s += w; } |
---|
606 | if (i) for (;;) { sum += *s; if (!(--i)) break; s += w; } |
---|
607 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
608 | } |
---|
609 | |
---|
610 | Real IdentityMatrix::trace() const |
---|
611 | { |
---|
612 | Real sum = *store * nrows_val; |
---|
613 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
614 | } |
---|
615 | |
---|
616 | |
---|
617 | Real BaseMatrix::trace() const |
---|
618 | { |
---|
619 | REPORT |
---|
620 | MatrixType Diag = MatrixType::Dg; Diag.SetDataLossOK(); |
---|
621 | GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(Diag); |
---|
622 | Real sum = gm->trace(); return sum; |
---|
623 | } |
---|
624 | |
---|
625 | void LogAndSign::operator*=(Real x) |
---|
626 | { |
---|
627 | if (x > 0.0) { log_val += log(x); } |
---|
628 | else if (x < 0.0) { log_val += log(-x); sign_val = -sign_val; } |
---|
629 | else sign_val = 0; |
---|
630 | } |
---|
631 | |
---|
632 | void LogAndSign::pow_eq(int k) |
---|
633 | { |
---|
634 | if (sign_val) |
---|
635 | { |
---|
636 | log_val *= k; |
---|
637 | if ( (k & 1) == 0 ) sign_val = 1; |
---|
638 | } |
---|
639 | } |
---|
640 | |
---|
641 | Real LogAndSign::value() const |
---|
642 | { |
---|
643 | Tracer et("LogAndSign::value"); |
---|
644 | if (log_val >= FloatingPointPrecision::LnMaximum()) |
---|
645 | Throw(OverflowException("Overflow in exponential")); |
---|
646 | return sign_val * exp(log_val); |
---|
647 | } |
---|
648 | |
---|
649 | LogAndSign::LogAndSign(Real f) |
---|
650 | { |
---|
651 | if (f == 0.0) { log_val = 0.0; sign_val = 0; return; } |
---|
652 | else if (f < 0.0) { sign_val = -1; f = -f; } |
---|
653 | else sign_val = 1; |
---|
654 | log_val = log(f); |
---|
655 | } |
---|
656 | |
---|
657 | LogAndSign DiagonalMatrix::log_determinant() const |
---|
658 | { |
---|
659 | REPORT |
---|
660 | int i = nrows_val; LogAndSign sum; Real* s = store; |
---|
661 | while (i--) sum *= *s++; |
---|
662 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
663 | } |
---|
664 | |
---|
665 | LogAndSign LowerTriangularMatrix::log_determinant() const |
---|
666 | { |
---|
667 | REPORT |
---|
668 | int i = nrows_val; LogAndSign sum; Real* s = store; int j = 2; |
---|
669 | // while (i--) { sum *= *s; s += j++; } |
---|
670 | if (i) for(;;) { sum *= *s; if (!(--i)) break; s += j++; } |
---|
671 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
672 | } |
---|
673 | |
---|
674 | LogAndSign UpperTriangularMatrix::log_determinant() const |
---|
675 | { |
---|
676 | REPORT |
---|
677 | int i = nrows_val; LogAndSign sum; Real* s = store; |
---|
678 | while (i) { sum *= *s; s += i--; } |
---|
679 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
680 | } |
---|
681 | |
---|
682 | LogAndSign IdentityMatrix::log_determinant() const |
---|
683 | { |
---|
684 | REPORT |
---|
685 | int i = nrows_val; LogAndSign sum; |
---|
686 | if (i > 0) { sum = *store; sum.PowEq(i); } |
---|
687 | ((GeneralMatrix&)*this).tDelete(); return sum; |
---|
688 | } |
---|
689 | |
---|
690 | LogAndSign BaseMatrix::log_determinant() const |
---|
691 | { |
---|
692 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
693 | LogAndSign sum = gm->log_determinant(); return sum; |
---|
694 | } |
---|
695 | |
---|
696 | LogAndSign GeneralMatrix::log_determinant() const |
---|
697 | { |
---|
698 | REPORT |
---|
699 | Tracer tr("log_determinant"); |
---|
700 | if (nrows_val != ncols_val) Throw(NotSquareException(*this)); |
---|
701 | CroutMatrix C(*this); return C.log_determinant(); |
---|
702 | } |
---|
703 | |
---|
704 | LogAndSign CroutMatrix::log_determinant() const |
---|
705 | { |
---|
706 | REPORT |
---|
707 | if (sing) return 0.0; |
---|
708 | int i = nrows_val; int dd = i+1; LogAndSign sum; Real* s = store; |
---|
709 | if (i) for(;;) |
---|
710 | { |
---|
711 | sum *= *s; |
---|
712 | if (!(--i)) break; |
---|
713 | s += dd; |
---|
714 | } |
---|
715 | if (!d) sum.ChangeSign(); return sum; |
---|
716 | |
---|
717 | } |
---|
718 | |
---|
719 | Real BaseMatrix::determinant() const |
---|
720 | { |
---|
721 | REPORT |
---|
722 | Tracer tr("determinant"); |
---|
723 | REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
724 | LogAndSign ld = gm->log_determinant(); |
---|
725 | return ld.Value(); |
---|
726 | } |
---|
727 | |
---|
728 | LinearEquationSolver::LinearEquationSolver(const BaseMatrix& bm) |
---|
729 | { |
---|
730 | gm = ( ((BaseMatrix&)bm).Evaluate() )->MakeSolver(); |
---|
731 | if (gm==&bm) { REPORT gm = gm->Image(); } |
---|
732 | // want a copy if *gm is actually bm |
---|
733 | else { REPORT gm->Protect(); } |
---|
734 | } |
---|
735 | |
---|
736 | ReturnMatrix BaseMatrix::sum_square_rows() const |
---|
737 | { |
---|
738 | REPORT |
---|
739 | GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
740 | int nr = gm->nrows(); |
---|
741 | ColumnVector ssq(nr); |
---|
742 | if (gm->size() == 0) { REPORT ssq = 0.0; } |
---|
743 | else |
---|
744 | { |
---|
745 | MatrixRow mr(gm, LoadOnEntry); |
---|
746 | for (int i = 1; i <= nr; ++i) |
---|
747 | { |
---|
748 | Real sum = 0.0; |
---|
749 | int s = mr.Storage(); |
---|
750 | Real* in = mr.Data(); |
---|
751 | while (s--) sum += square(*in++); |
---|
752 | ssq(i) = sum; |
---|
753 | mr.Next(); |
---|
754 | } |
---|
755 | } |
---|
756 | gm->tDelete(); |
---|
757 | ssq.release(); return ssq.for_return(); |
---|
758 | } |
---|
759 | |
---|
760 | ReturnMatrix BaseMatrix::sum_square_columns() const |
---|
761 | { |
---|
762 | REPORT |
---|
763 | GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
764 | int nr = gm->nrows(); int nc = gm->ncols(); |
---|
765 | RowVector ssq(nc); ssq = 0.0; |
---|
766 | if (gm->size() != 0) |
---|
767 | { |
---|
768 | MatrixRow mr(gm, LoadOnEntry); |
---|
769 | for (int i = 1; i <= nr; ++i) |
---|
770 | { |
---|
771 | int s = mr.Storage(); |
---|
772 | Real* in = mr.Data(); Real* out = ssq.data() + mr.Skip(); |
---|
773 | while (s--) *out++ += square(*in++); |
---|
774 | mr.Next(); |
---|
775 | } |
---|
776 | } |
---|
777 | gm->tDelete(); |
---|
778 | ssq.release(); return ssq.for_return(); |
---|
779 | } |
---|
780 | |
---|
781 | ReturnMatrix BaseMatrix::sum_rows() const |
---|
782 | { |
---|
783 | REPORT |
---|
784 | GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
785 | int nr = gm->nrows(); |
---|
786 | ColumnVector sum_vec(nr); |
---|
787 | if (gm->size() == 0) { REPORT sum_vec = 0.0; } |
---|
788 | else |
---|
789 | { |
---|
790 | MatrixRow mr(gm, LoadOnEntry); |
---|
791 | for (int i = 1; i <= nr; ++i) |
---|
792 | { |
---|
793 | Real sum = 0.0; |
---|
794 | int s = mr.Storage(); |
---|
795 | Real* in = mr.Data(); |
---|
796 | while (s--) sum += *in++; |
---|
797 | sum_vec(i) = sum; |
---|
798 | mr.Next(); |
---|
799 | } |
---|
800 | } |
---|
801 | gm->tDelete(); |
---|
802 | sum_vec.release(); return sum_vec.for_return(); |
---|
803 | } |
---|
804 | |
---|
805 | ReturnMatrix BaseMatrix::sum_columns() const |
---|
806 | { |
---|
807 | REPORT |
---|
808 | GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
809 | int nr = gm->nrows(); int nc = gm->ncols(); |
---|
810 | RowVector sum_vec(nc); sum_vec = 0.0; |
---|
811 | if (gm->size() != 0) |
---|
812 | { |
---|
813 | MatrixRow mr(gm, LoadOnEntry); |
---|
814 | for (int i = 1; i <= nr; ++i) |
---|
815 | { |
---|
816 | int s = mr.Storage(); |
---|
817 | Real* in = mr.Data(); Real* out = sum_vec.data() + mr.Skip(); |
---|
818 | while (s--) *out++ += *in++; |
---|
819 | mr.Next(); |
---|
820 | } |
---|
821 | } |
---|
822 | gm->tDelete(); |
---|
823 | sum_vec.release(); return sum_vec.for_return(); |
---|
824 | } |
---|
825 | |
---|
826 | |
---|
827 | #ifdef use_namespace |
---|
828 | } |
---|
829 | #endif |
---|
830 | |
---|
831 | |
---|
832 | ///} |
---|