source: ntrip/trunk/BNC/newmat/newmat8.cpp@ 10113

Last change on this file since 10113 was 9434, checked in by stuerze, 3 years ago

newmat is reverted to its stable version 10 because version 11beta was not working with mac compilers
File size: 21.9 KB
Line 
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
18namespace 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
31void 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
103void 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
154inline Real square(Real x) { return x*x; }
155
156Real 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
164Real 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
172Real 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
212static void NullMatrixError(const GeneralMatrix* gm)
213{
214 ((GeneralMatrix&)*gm).tDelete();
215 Throw(ProgramException("Maximum or minimum of null matrix"));
216}
217
218Real 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
227Real 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
238Real 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
247Real 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
258Real 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
267Real 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
277Real 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
286Real 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
296Real 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
311Real 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
326Real 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
341Real 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
356Real 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
364Real 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
372Real 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
380Real 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
388Real 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
401Real 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
414Real IdentityMatrix::sum_absolute_value() const
415 { REPORT return fabs(trace()); } // no need to do tDelete?
416
417Real 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
430Real IdentityMatrix::sum_square() const
431{
432 Real sum = *store * *store * nrows_val;
433 ((GeneralMatrix&)*this).tDelete(); return sum;
434}
435
436
437Real BaseMatrix::sum_square() const
438{
439 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
440 Real s = gm->sum_square(); return s;
441}
442
443Real BaseMatrix::norm_Frobenius() const
444 { REPORT return sqrt(sum_square()); }
445
446Real BaseMatrix::sum_absolute_value() const
447{
448 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
449 Real s = gm->sum_absolute_value(); return s;
450}
451
452Real BaseMatrix::sum() const
453{
454 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
455 Real s = gm->sum(); return s;
456}
457
458Real BaseMatrix::maximum_absolute_value() const
459{
460 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
461 Real s = gm->maximum_absolute_value(); return s;
462}
463
464Real 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
470Real 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
476Real BaseMatrix::minimum_absolute_value() const
477{
478 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
479 Real s = gm->minimum_absolute_value(); return s;
480}
481
482Real 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
488Real 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
494Real BaseMatrix::maximum() const
495{
496 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
497 Real s = gm->maximum(); return s;
498}
499
500Real BaseMatrix::maximum1(int& i) const
501{
502 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
503 Real s = gm->maximum1(i); return s;
504}
505
506Real 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
512Real BaseMatrix::minimum() const
513{
514 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
515 Real s = gm->minimum(); return s;
516}
517
518Real BaseMatrix::minimum1(int& i) const
519{
520 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
521 Real s = gm->minimum1(i); return s;
522}
523
524Real 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
530Real 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
544Real 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
556Real 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
564Real 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
573Real 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
582Real 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
590Real 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
600Real 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
610Real IdentityMatrix::trace() const
611{
612 Real sum = *store * nrows_val;
613 ((GeneralMatrix&)*this).tDelete(); return sum;
614}
615
616
617Real 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
625void 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
632void 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
641Real 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
649LogAndSign::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
657LogAndSign 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
665LogAndSign 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
674LogAndSign 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
682LogAndSign 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
690LogAndSign BaseMatrix::log_determinant() const
691{
692 REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
693 LogAndSign sum = gm->log_determinant(); return sum;
694}
695
696LogAndSign 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
704LogAndSign 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
719Real 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
728LinearEquationSolver::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
736ReturnMatrix 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
760ReturnMatrix 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
781ReturnMatrix 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
805ReturnMatrix 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///}
Note: See TracBrowser for help on using the repository browser.