newmat8.cpp
Go to the documentation of this file.
1 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
297 {
298  REPORT
299  if (storage == 0) NullMatrixError(this);
300  Real maxval = 0.0; int nr = Nrows();
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 
312 {
313  REPORT
314  if (storage == 0) NullMatrixError(this);
315  Real minval = FloatingPointPrecision::Maximum(); int nr = Nrows();
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();
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();
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 
357 {
358  REPORT
360  i = k / Ncols(); j = k - i * Ncols(); i++; j++;
361  return m;
362 }
363 
365 {
366  REPORT
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 
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 
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 
415  { REPORT return fabs(trace()); } // no need to do tDelete?
416 
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 
431 {
432  Real sum = *store * *store * nrows_val;
433  ((GeneralMatrix&)*this).tDelete(); return sum;
434 }
435 
436 
438 {
439  REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
440  Real s = gm->sum_square(); return s;
441 }
442 
444  { REPORT return sqrt(sum_square()); }
445 
447 {
448  REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
449  Real s = gm->sum_absolute_value(); return s;
450 }
451 
453 {
454  REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
455  Real s = gm->sum(); return s;
456 }
457 
459 {
460  REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
461  Real s = gm->maximum_absolute_value(); return s;
462 }
463 
465 {
466  REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
467  Real s = gm->maximum_absolute_value1(i); return s;
468 }
469 
471 {
472  REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
473  Real s = gm->maximum_absolute_value2(i, j); return s;
474 }
475 
477 {
478  REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
479  Real s = gm->minimum_absolute_value(); return s;
480 }
481 
483 {
484  REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
485  Real s = gm->minimum_absolute_value1(i); return s;
486 }
487 
489 {
490  REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
491  Real s = gm->minimum_absolute_value2(i, j); return s;
492 }
493 
495 {
496  REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
497  Real s = gm->maximum(); return s;
498 }
499 
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 
513 {
514  REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
515  Real s = gm->minimum(); return s;
516 }
517 
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");
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 
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 
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 
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 
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 
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 
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 
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 
611 {
612  Real sum = *store * nrows_val;
613  ((GeneralMatrix&)*this).tDelete(); return sum;
614 }
615 
616 
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 
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 
633 {
634  if (sign_val)
635  {
636  log_val *= k;
637  if ( (k & 1) == 0 ) sign_val = 1;
638  }
639 }
640 
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 
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 
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 
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 
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 
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 
691 {
692  REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
693  LogAndSign sum = gm->log_determinant(); return sum;
694 }
695 
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 
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 
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 
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 
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 
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 
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 
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 
virtual Real minimum_absolute_value() const
Definition: newmat8.cpp:476
Real minimum2(int &i, int &j) const
Definition: newmat8.cpp:341
void ludcmp()
Definition: newmat8.cpp:31
LogAndSign log_determinant() const
Definition: newmat8.cpp:682
Real norm_Frobenius() const
Definition: newmat8.cpp:443
int Storage()
Definition: newmatrc.h:101
Real trace() const
Definition: newmat8.cpp:600
Real minimum2(int &i, int &j) const
Definition: newmat8.cpp:380
Miscellaneous exception (details in character string).
Definition: newmat.h:1947
static Real Maximum()
Definition: precisio.h:235
Real maximum_absolute_value2(int &i, int &j) const
Definition: newmat8.cpp:296
Real Value() const
Definition: newmat.h:64
virtual Real sum_square() const
Definition: newmat8.cpp:437
Real trace() const
Definition: newmat8.cpp:590
ReturnMatrix for_return() const
Definition: newmat4.cpp:249
virtual Real maximum() const
Definition: newmat8.cpp:494
LogAndSign log_determinant() const
Definition: newmat8.cpp:704
Real trace() const
Definition: newmat8.cpp:582
static void NullMatrixError(const GeneralMatrix *gm)
Definition: newmat8.cpp:212
Real square(Real x)
Definition: newmat8.cpp:154
void operator*=(Real)
multiply by a real
Definition: newmat8.cpp:625
static Real LnMaximum()
Definition: precisio.h:238
Real minimum_absolute_value1(int &i) const
Definition: newmat8.cpp:247
virtual Real maximum2(int &i, int &j) const
Definition: newmat8.cpp:506
void pow_eq(int k)
raise to power of k
Definition: newmat8.cpp:632
virtual Real trace() const
Definition: newmat8.cpp:617
virtual Real minimum2(int &i, int &j) const
Definition: newmat8.cpp:524
double Real
Definition: include.h:307
ReturnMatrix sum_rows() const
Definition: newmat8.cpp:781
Real minimum_absolute_value2(int &i, int &j) const
Definition: newmat8.cpp:364
void PowEq(int k)
Definition: newmat.h:56
GeneralMatrix * Evaluate(MatrixType mt=MatrixTypeUnSp)
Definition: newmat5.cpp:85
int Skip()
Definition: newmatrc.h:100
Real sum() const
Definition: newmat8.cpp:172
virtual Real minimum_absolute_value2(int &i, int &j) const
Definition: newmat8.cpp:488
LogAndSign()
Definition: newmat.h:52
#define REPORT
Definition: newmat8.cpp:25
ReturnMatrix sum_columns() const
Definition: newmat8.cpp:805
virtual Real maximum_absolute_value1(int &i) const
Definition: newmat8.cpp:464
Real sum_square() const
Definition: newmat8.cpp:388
Real maximum_absolute_value() const
Definition: newmat8.cpp:218
Real maximum() const
Definition: newmat8.cpp:258
Real * data()
Definition: newmat.h:502
virtual Real minimum_absolute_value1(int &i) const
Definition: newmat8.cpp:482
Real trace() const
Definition: newmat8.cpp:564
A matrix is not square exception.
Definition: newmat.h:1981
virtual Real minimum() const
Definition: newmat8.cpp:512
LinearEquationSolver(const BaseMatrix &bm)
Definition: newmat8.cpp:728
Real minimum_absolute_value() const
Definition: newmat8.cpp:238
Real sum_absolute_value() const
Definition: newmat8.cpp:401
virtual Real minimum1(int &i) const
Definition: newmat8.cpp:518
FloatVector * d
LogAndSign log_determinant() const
Definition: newmat8.cpp:696
Real maximum1(int &i) const
Definition: newmat8.cpp:267
Real value() const
the value
Definition: newmat8.cpp:641
Real sum_absolute_value() const
Definition: newmat8.cpp:164
Real maximum2(int &i, int &j) const
Definition: newmat8.cpp:326
Real * Data()
Definition: newmatrc.h:99
Real trace() const
Definition: newmat8.cpp:610
void lubksb(Real *, int=0)
Definition: newmat8.cpp:103
Real Maximum1(Real r, int &i)
Definition: newmat2.cpp:616
virtual Real sum_absolute_value() const
Definition: newmat8.cpp:446
Real trace() const
Definition: newmat8.cpp:544
Real sum_square() const
Definition: newmat8.cpp:156
Real maximum_absolute_value1(int &i) const
Definition: newmat8.cpp:227
virtual Real maximum_absolute_value() const
Definition: newmat8.cpp:458
Singular matrix exception.
Definition: newmat.h:1931
Real maximum2(int &i, int &j) const
Definition: newmat8.cpp:372
Real dotproduct(const Matrix &A, const Matrix &B)
Definition: newmat8.cpp:530
Real maximum_absolute_value2(int &i, int &j) const
Definition: newmat8.cpp:356
Base of the matrix classes.
Definition: newmat.h:292
Real sum_square(const BaseMatrix &B)
Definition: newmat.h:2094
LogAndSign log_determinant() const
Definition: newmat8.cpp:665
#define Throw(E)
Definition: myexcept.h:191
ReturnMatrix sum_square_rows() const
Definition: newmat8.cpp:736
LogAndSign log_determinant() const
Definition: newmat8.cpp:657
Real sum_square() const
Definition: newmat8.cpp:430
The usual rectangular matrix.
Definition: newmat.h:625
Real * store
Definition: newmat.h:454
Real MaximumAbsoluteValue1(Real r, int &i)
Definition: newmat2.cpp:596
Incompatible dimensions exception.
Definition: newmat.h:1998
virtual Real maximum_absolute_value2(int &i, int &j) const
Definition: newmat8.cpp:470
Real determinant() const
Definition: newmat8.cpp:719
FloatVector FloatVector * a
LogAndSign log_determinant() const
Definition: newmat8.cpp:674
Real minimum_absolute_value2(int &i, int &j) const
Definition: newmat8.cpp:311
Real trace(const BaseMatrix &B)
Definition: newmat.h:2099
Real trace() const
Definition: newmat8.cpp:556
Real trace() const
Definition: newmat8.cpp:573
Real overflow exception.
Definition: newmat.h:1939
ReturnMatrix sum_square_columns() const
Definition: newmat8.cpp:760
Real minimum() const
Definition: newmat8.cpp:277
Real sum(const BaseMatrix &B)
Definition: newmat.h:2105
int nrows() const
Definition: newmat.h:499
void ChangeSign()
Definition: newmat.h:57
Real sum() const
Definition: newmat8.cpp:417
Row vector.
Definition: newmat.h:953
Real MinimumAbsoluteValue1(Real r, int &i)
Definition: newmat2.cpp:606
int storage
Definition: newmat.h:453
Column vector.
Definition: newmat.h:1008
The classes for matrices that can contain data are derived from this.
Definition: newmat.h:447
Real sum_absolute_value() const
Definition: newmat8.cpp:414
virtual Real sum() const
Definition: newmat8.cpp:452
void Next()
Definition: newmatrc.h:161
void release()
Definition: newmat.h:517
Real minimum1(int &i) const
Definition: newmat8.cpp:286
Real Minimum1(Real r, int &i)
Definition: newmat2.cpp:626
virtual Real maximum1(int &i) const
Definition: newmat8.cpp:500
virtual LogAndSign log_determinant() const
Definition: newmat8.cpp:690


kni
Author(s): Martin Günther
autogenerated on Fri Jan 3 2020 04:01:16