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
ludcmp()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
lubksb(Real * B,int mini)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
square(Real x)154 inline Real square(Real x) { return x*x; }
155
sum_square() const156 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
sum_absolute_value() const164 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
sum() const172 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
NullMatrixError(const GeneralMatrix * gm)212 static void NullMatrixError(const GeneralMatrix* gm)
213 {
214 ((GeneralMatrix&)*gm).tDelete();
215 Throw(ProgramException("Maximum or minimum of null matrix"));
216 }
217
maximum_absolute_value() const218 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
maximum_absolute_value1(int & i) const227 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
minimum_absolute_value() const238 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
minimum_absolute_value1(int & i) const247 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
maximum() const258 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
maximum1(int & i) const267 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
minimum() const277 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
minimum1(int & i) const286 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
maximum_absolute_value2(int & i,int & j) const296 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
minimum_absolute_value2(int & i,int & j) const311 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
maximum2(int & i,int & j) const326 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
minimum2(int & i,int & j) const341 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
maximum_absolute_value2(int & i,int & j) const356 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
minimum_absolute_value2(int & i,int & j) const364 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
maximum2(int & i,int & j) const372 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
minimum2(int & i,int & j) const380 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
sum_square() const388 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
sum_absolute_value() const401 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
sum_absolute_value() const414 Real IdentityMatrix::sum_absolute_value() const
415 { REPORT return fabs(trace()); } // no need to do tDelete?
416
sum() const417 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
sum_square() const430 Real IdentityMatrix::sum_square() const
431 {
432 Real sum = *store * *store * nrows_val;
433 ((GeneralMatrix&)*this).tDelete(); return sum;
434 }
435
436
sum_square() const437 Real BaseMatrix::sum_square() const
438 {
439 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
440 Real s = gm->sum_square(); return s;
441 }
442
norm_Frobenius() const443 Real BaseMatrix::norm_Frobenius() const
444 { REPORT return sqrt(sum_square()); }
445
sum_absolute_value() const446 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
sum() const452 Real BaseMatrix::sum() const
453 {
454 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
455 Real s = gm->sum(); return s;
456 }
457
maximum_absolute_value() const458 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
maximum_absolute_value1(int & i) const464 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
maximum_absolute_value2(int & i,int & j) const470 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
minimum_absolute_value() const476 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
minimum_absolute_value1(int & i) const482 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
minimum_absolute_value2(int & i,int & j) const488 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
maximum() const494 Real BaseMatrix::maximum() const
495 {
496 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
497 Real s = gm->maximum(); return s;
498 }
499
maximum1(int & i) const500 Real BaseMatrix::maximum1(int& i) const
501 {
502 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
503 Real s = gm->maximum1(i); return s;
504 }
505
maximum2(int & i,int & j) const506 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
minimum() const512 Real BaseMatrix::minimum() const
513 {
514 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
515 Real s = gm->minimum(); return s;
516 }
517
minimum1(int & i) const518 Real BaseMatrix::minimum1(int& i) const
519 {
520 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
521 Real s = gm->minimum1(i); return s;
522 }
523
minimum2(int & i,int & j) const524 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
dotproduct(const Matrix & A,const Matrix & B)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
trace() const544 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
trace() const556 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
trace() const564 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
trace() const573 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
trace() const582 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
trace() const590 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
trace() const600 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
trace() const610 Real IdentityMatrix::trace() const
611 {
612 Real sum = *store * nrows_val;
613 ((GeneralMatrix&)*this).tDelete(); return sum;
614 }
615
616
trace() const617 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
operator *=(Real x)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
pow_eq(int k)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
value() const641 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
LogAndSign(Real f)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
log_determinant() const657 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
log_determinant() const665 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
log_determinant() const674 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
log_determinant() const682 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
log_determinant() const690 LogAndSign BaseMatrix::log_determinant() const
691 {
692 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
693 LogAndSign sum = gm->log_determinant(); return sum;
694 }
695
log_determinant() const696 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
log_determinant() const704 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
determinant() const719 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
LinearEquationSolver(const BaseMatrix & bm)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
sum_square_rows() const736 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
sum_square_columns() const760 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
sum_rows() const781 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
sum_columns() const805 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 ///}
833