Logo ROOT   6.21/01
Reference Guide
TMatrixTBase.cxx
Go to the documentation of this file.
1 // @(#)root/matrix:$Id: 2d00df45ce4c38c7ea0930d6b520cbf4cfb9152e $
2 // Authors: Fons Rademakers, Eddy Offermann Nov 2003
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /** \class TMatrixTBase
13  \ingroup Matrix
14  \brief Linear Algebra Package
15 
16  ### Linear Algebra Package
17 
18  The present package implements all the basic algorithms dealing
19  with vectors, matrices, matrix columns, rows, diagonals, etc.
20  In addition eigen-Vector analysis and several matrix decomposition
21  have been added (LU,QRH,Cholesky,Bunch-Kaufman and SVD) .
22  The decompositions are used in matrix inversion, equation solving.
23 
24  For a dense matrix, elements are arranged in memory in a ROW-wise
25  fashion . For (n x m) matrices where n*m <=kSizeMax (=25 currently)
26  storage space is available on the stack, thus avoiding expensive
27  allocation/deallocation of heap space . However, this introduces of
28  course kSizeMax overhead for each matrix object . If this is an
29  issue recompile with a new appropriate value (>=0) for kSizeMax
30 
31  Sparse matrices are also stored in row-wise fashion but additional
32  row/column information is stored, see TMatrixTSparse source for
33  additional details .
34 
35  Another way to assign and store matrix data is through Use
36  see for instance stressLinear.cxx file .
37 
38  Unless otherwise specified, matrix and vector indices always start
39  with 0, spanning up to the specified limit-1. However, there are
40  constructors to which one can specify aribtrary lower and upper
41  bounds, e.g. TMatrixD m(1,10,1,5) defines a matrix that ranges
42  from 1..10, 1..5 (a(1,1)..a(10,5)).
43 
44  The present package provides all facilities to completely AVOID
45  returning matrices. Use "TMatrixD A(TMatrixD::kTransposed,B);"
46  and other fancy constructors as much as possible. If one really needs
47  to return a matrix, return a TMatrixTLazy object instead. The
48  conversion is completely transparent to the end user, e.g.
49  "TMatrixT m = THaarMatrixT(5);" and _is_ efficient.
50 
51  Since TMatrixT et al. are fully integrated in ROOT, they of course
52  can be stored in a ROOT database.
53 
54  For usage examples see $ROOTSYS/test/stressLinear.cxx
55 
56  ### Acknowledgements
57 
58  1. Oleg E. Kiselyov
59  First implementations were based on the his code . We have diverged
60  quite a bit since then but the ideas/code for lazy matrix and
61  "nested function" are 100% his .
62  You can see him and his code in action at http://okmij.org/ftp
63  2. Chris R. Birchenhall,
64  We adapted his idea of the implementation for the decomposition
65  classes instead of our messy installation of matrix inversion
66  His installation of matrix condition number, using an iterative
67  scheme using the Hage algorithm is worth looking at !
68  Chris has a nice writeup (matdoc.ps) on his matrix classes at
69  ftp://ftp.mcc.ac.uk/pub/matclass/
70  3. Mark Fischler and Steven Haywood of CLHEP
71  They did the slave labor of spelling out all sub-determinants
72  for Cramer inversion of (4x4),(5x5) and (6x6) matrices
73  The stack storage for small matrices was also taken from them
74  4. Roldan Pozo of TNT (http://math.nist.gov/tnt/)
75  He converted the EISPACK routines for the eigen-vector analysis to
76  C++ . We started with his implementation
77  5. Siegmund Brandt (http://siux00.physik.uni-siegen.de/~brandt/datan
78  We adapted his (very-well) documented SVD routines
79 
80 ### How to efficiently use this package
81 
82 #### 1. Never return complex objects (matrices or vectors)
83  Danger: For example, when the following snippet:
84 ~~~
85  TMatrixD foo(int n)
86  {
87  TMatrixD foom(n,n); fill_in(foom); return foom;
88  }
89  TMatrixD m = foo(5);
90 ~~~
91  runs, it constructs matrix foo:foom, copies it onto stack as a
92  return value and destroys foo:foom. Return value (a matrix)
93  from foo() is then copied over to m (via a copy constructor),
94  and the return value is destroyed. So, the matrix constructor is
95  called 3 times and the destructor 2 times. For big matrices,
96  the cost of multiple constructing/copying/destroying of objects
97  may be very large. *Some* optimized compilers can cut down on 1
98  copying/destroying, but still it leaves at least two calls to
99  the constructor. Note, TMatrixDLazy (see below) can construct
100  TMatrixD m "inplace", with only a _single_ call to the
101  constructor.
102 
103 #### 2. Use "two-address instructions"
104 ~~~
105  "void TMatrixD::operator += (const TMatrixD &B);"
106 ~~~
107  as much as possible.
108  That is, to add two matrices, it's much more efficient to write
109 ~~~
110  A += B;
111 ~~~
112  than
113 ~~~
114  TMatrixD C = A + B;
115 ~~~
116  (if both operand should be preserved,
117  TMatrixD C = A; C += B;
118  is still better).
119 
120 #### 3. Use glorified constructors when returning of an object seems inevitable:
121 ~~~
122  "TMatrixD A(TMatrixD::kTransposed,B);"
123  "TMatrixD C(A,TMatrixD::kTransposeMult,B);"
124 ~~~
125 
126  like in the following snippet (from $ROOTSYS/test/vmatrix.cxx)
127  that verifies that for an orthogonal matrix T, T'T = TT' = E.
128 
129 ~~~
130  TMatrixD haar = THaarMatrixD(5);
131  TMatrixD unit(TMatrixD::kUnit,haar);
132  TMatrixD haar_t(TMatrixD::kTransposed,haar);
133  TMatrixD hth(haar,TMatrixD::kTransposeMult,haar);
134  TMatrixD hht(haar,TMatrixD::kMult,haar_t);
135  TMatrixD hht1 = haar; hht1 *= haar_t;
136  VerifyMatrixIdentity(unit,hth);
137  VerifyMatrixIdentity(unit,hht);
138  VerifyMatrixIdentity(unit,hht1);
139 ~~~
140 
141 #### 4. Accessing row/col/diagonal of a matrix without much fuss
142  (and without moving a lot of stuff around):
143 
144 ~~~
145  TMatrixD m(n,n); TVectorD v(n); TMatrixDDiag(m) += 4;
146  v = TMatrixDRow(m,0);
147  TMatrixDColumn m1(m,1); m1(2) = 3; // the same as m(2,1)=3;
148 ~~~
149  Note, constructing of, say, TMatrixDDiag does *not* involve any
150  copying of any elements of the source matrix.
151 
152 #### 5. It's possible (and encouraged) to use "nested" functions
153  For example, creating of a Hilbert matrix can be done as follows:
154 
155 ~~~
156  void foo(const TMatrixD &m)
157  {
158  TMatrixD m1(TMatrixD::kZero,m);
159  struct MakeHilbert : public TElementPosActionD {
160  void Operation(Double_t &element)
161  { element = 1./(fI+fJ-1); }
162  };
163  m1.Apply(MakeHilbert());
164  }
165 ~~~
166 
167  of course, using a special method THilbertMatrixD() is
168  still more optimal, but not by a whole lot. And that's right,
169  class MakeHilbert is declared *within* a function and local to
170  that function. It means one can define another MakeHilbert class
171  (within another function or outside of any function, that is, in
172  the global scope), and it still will be OK. Note, this currently
173  is not yet supported by the interpreter CINT.
174 
175  Another example is applying of a simple function to each matrix element:
176 
177 ~~~
178  void foo(TMatrixD &m,TMatrixD &m1)
179  {
180  typedef double (*dfunc_t)(double);
181  class ApplyFunction : public TElementActionD {
182  dfunc_t fFunc;
183  void Operation(Double_t &element)
184  { element=fFunc(element); }
185  public:
186  ApplyFunction(dfunc_t func):fFunc(func) {}
187  };
188  ApplyFunction x(TMath::Sin);
189  m.Apply(x);
190  }
191 ~~~
192 
193  Validation code $ROOTSYS/test/vmatrix.cxx and vvector.cxx contain
194  a few more examples of that kind.
195 
196 #### 6. Lazy matrices:
197  instead of returning an object return a "recipe"
198  how to make it. The full matrix would be rolled out only when
199  and where it's needed:
200 ~~~
201  TMatrixD haar = THaarMatrixD(5);
202 ~~~
203  THaarMatrixD() is a *class*, not a simple function. However
204  similar this looks to a returning of an object (see note #1
205  above), it's dramatically different. THaarMatrixD() constructs a
206  TMatrixDLazy, an object of just a few bytes long. A special
207  "TMatrixD(const TMatrixDLazy &recipe)" constructor follows the
208  recipe and makes the matrix haar() right in place. No matrix
209  element is moved whatsoever!
210 
211 */
212 
213 ////////////////////////////////////////////////////////////////////////////////
214 ///
215 /// TMatrixTBase
216 ///
217 /// Template of base class in the linear algebra package
218 ///
219 /// matrix properties are stored here, however the data storage is part
220 /// of the derived classes
221 
222 #include "TMatrixTBase.h"
223 #include "TVectorT.h"
224 #include "TROOT.h"
225 #include "TClass.h"
226 #include "TMath.h"
227 #include <limits.h>
228 
230 
232 
233 
234 ////////////////////////////////////////////////////////////////////////////////
235 /// Lexical sort on array data using indices first and second
236 
237 template<class Element>
239 {
240  const int incs[] = {1,5,19,41,109,209,505,929,2161,3905,8929,16001,INT_MAX};
241 
242  Int_t kinc = 0;
243  while (incs[kinc] <= n/2)
244  kinc++;
245  kinc -= 1;
246 
247  // incs[kinc] is the greatest value in the sequence that is also <= n/2.
248  // If n == {0,1}, kinc == -1 and so no sort will take place.
249 
250  for( ; kinc >= 0; kinc--) {
251  const Int_t inc = incs[kinc];
252 
253  for (Int_t k = inc; k < n; k++) {
254  const Element tmp = data[k];
255  const Int_t fi = first [k];
256  const Int_t se = second[k];
257  Int_t j;
258  for (j = k; j >= inc; j -= inc) {
259  if ( fi < first[j-inc] || (fi == first[j-inc] && se < second[j-inc]) ) {
260  data [j] = data [j-inc];
261  first [j] = first [j-inc];
262  second[j] = second[j-inc];
263  } else
264  break;
265  }
266  data [j] = tmp;
267  first [j] = fi;
268  second[j] = se;
269  }
270  }
271 }
272 
273 ////////////////////////////////////////////////////////////////////////////////
274 /// Lexical sort on array data using indices first and second
275 
276 template<class Element>
278  Int_t *second,Int_t swapSecond,Int_t *index)
279 {
280  const int incs[] = {1,5,19,41,109,209,505,929,2161,3905,8929,16001,INT_MAX};
281 
282  Int_t kinc = 0;
283  while (incs[kinc] <= n/2)
284  kinc++;
285  kinc -= 1;
286 
287  // incs[kinc] is the greatest value in the sequence that is also less
288  // than n/2.
289 
290  for( ; kinc >= 0; kinc--) {
291  const Int_t inc = incs[kinc];
292 
293  if ( !swapFirst && !swapSecond ) {
294  for (Int_t k = inc; k < n; k++) {
295  // loop over all subarrays defined by the current increment
296  const Int_t ktemp = index[k];
297  const Int_t fi = first [ktemp];
298  const Int_t se = second[ktemp];
299  // Insert element k into the sorted subarray
300  Int_t j;
301  for (j = k; j >= inc; j -= inc) {
302  // Loop over the elements in the current subarray
303  if (fi < first[index[j-inc]] || (fi == first[index[j-inc]] && se < second[index[j-inc]])) {
304  // Swap elements j and j - inc, implicitly use the fact
305  // that ktemp hold element j to avoid having to assign to
306  // element j-inc
307  index[j] = index[j-inc];
308  } else {
309  // There are no more elements in this sorted subarray which
310  // are less than element j
311  break;
312  }
313  } // End loop over the elements in the current subarray
314  // Move index[j] out of temporary storage
315  index[j] = ktemp;
316  // The element has been inserted into the subarray.
317  } // End loop over all subarrays defined by the current increment
318  } else if ( swapSecond && !swapFirst ) {
319  for (Int_t k = inc; k < n; k++) {
320  const Int_t ktemp = index[k];
321  const Int_t fi = first [ktemp];
322  const Int_t se = second[k];
323  Int_t j;
324  for (j = k; j >= inc; j -= inc) {
325  if (fi < first[index[j-inc]] || (fi == first[index[j-inc]] && se < second[j-inc])) {
326  index [j] = index[j-inc];
327  second[j] = second[j-inc];
328  } else {
329  break;
330  }
331  }
332  index[j] = ktemp;
333  second[j] = se;
334  }
335  } else if (swapFirst && !swapSecond) {
336  for (Int_t k = inc; k < n; k++ ) {
337  const Int_t ktemp = index[k];
338  const Int_t fi = first[k];
339  const Int_t se = second[ktemp];
340  Int_t j;
341  for (j = k; j >= inc; j -= inc) {
342  if ( fi < first[j-inc] || (fi == first[j-inc] && se < second[ index[j-inc]])) {
343  index[j] = index[j-inc];
344  first[j] = first[j-inc];
345  } else {
346  break;
347  }
348  }
349  index[j] = ktemp;
350  first[j] = fi;
351  }
352  } else { // Swap both
353  for (Int_t k = inc; k < n; k++ ) {
354  const Int_t ktemp = index[k];
355  const Int_t fi = first [k];
356  const Int_t se = second[k];
357  Int_t j;
358  for (j = k; j >= inc; j -= inc) {
359  if ( fi < first[j-inc] || (fi == first[j-inc] && se < second[j-inc])) {
360  index [j] = index [j-inc];
361  first [j] = first [j-inc];
362  second[j] = second[j-inc];
363  } else {
364  break;
365  }
366  }
367  index[j] = ktemp;
368  first[j] = fi;
369  second[j] = se;
370  }
371  }
372  }
373 }
374 
375 ////////////////////////////////////////////////////////////////////////////////
376 /// Copy array data to matrix . It is assumed that array is of size >= fNelems
377 /// (=)))) fNrows*fNcols
378 /// option indicates how the data is stored in the array:
379 /// option =
380 /// - 'F' : column major (Fortran) m[i][j] = array[i+j*fNrows]
381 /// - else : row major (C) m[i][j] = array[i*fNcols+j] (default)
382 
383 template<class Element>
385 {
386  R__ASSERT(IsValid());
387 
388  TString opt = option;
389  opt.ToUpper();
390 
391  Element *elem = GetMatrixArray();
392  if (opt.Contains("F")) {
393  for (Int_t irow = 0; irow < fNrows; irow++) {
394  const Int_t off1 = irow*fNcols;
395  Int_t off2 = 0;
396  for (Int_t icol = 0; icol < fNcols; icol++) {
397  elem[off1+icol] = data[off2+irow];
398  off2 += fNrows;
399  }
400  }
401  }
402  else
403  memcpy(elem,data,fNelems*sizeof(Element));
404 
405  return *this;
406 }
407 
408 ////////////////////////////////////////////////////////////////////////////////
409 /// Check whether matrix is symmetric
410 
411 template<class Element>
413 {
414  R__ASSERT(IsValid());
415 
416  if ((fNrows != fNcols) || (fRowLwb != fColLwb))
417  return kFALSE;
418 
419  const Element * const elem = GetMatrixArray();
420  for (Int_t irow = 0; irow < fNrows; irow++) {
421  const Int_t rowOff = irow*fNcols;
422  Int_t colOff = 0;
423  for (Int_t icol = 0; icol < irow; icol++) {
424  if (elem[rowOff+icol] != elem[colOff+irow])
425  return kFALSE;
426  colOff += fNrows;
427  }
428  }
429  return kTRUE;
430 }
431 
432 ////////////////////////////////////////////////////////////////////////////////
433 /// Copy matrix data to array . It is assumed that array is of size >= fNelems
434 /// (=)))) fNrows*fNcols
435 /// option indicates how the data is stored in the array:
436 /// option =
437 /// - 'F' : column major (Fortran) array[i+j*fNrows] = m[i][j]
438 /// - else : row major (C) array[i*fNcols+j] = m[i][j] (default)
439 
440 template<class Element>
441 void TMatrixTBase<Element>::GetMatrix2Array(Element *data,Option_t *option) const
442 {
443  R__ASSERT(IsValid());
444 
445  TString opt = option;
446  opt.ToUpper();
447 
448  const Element * const elem = GetMatrixArray();
449  if (opt.Contains("F")) {
450  for (Int_t irow = 0; irow < fNrows; irow++) {
451  const Int_t off1 = irow*fNcols;
452  Int_t off2 = 0;
453  for (Int_t icol = 0; icol < fNcols; icol++) {
454  data[off2+irow] = elem[off1+icol];
455  off2 += fNrows;
456  }
457  }
458  }
459  else
460  memcpy(data,elem,fNelems*sizeof(Element));
461 }
462 
463 ////////////////////////////////////////////////////////////////////////////////
464 /// Copy n elements from array v to row rown starting at column coln
465 
466 template<class Element>
468 {
469  const Int_t arown = rown-fRowLwb;
470  const Int_t acoln = coln-fColLwb;
471  const Int_t nr = (n > 0) ? n : fNcols;
472 
473  if (gMatrixCheck) {
474  if (arown >= fNrows || arown < 0) {
475  Error("InsertRow","row %d out of matrix range",rown);
476  return *this;
477  }
478 
479  if (acoln >= fNcols || acoln < 0) {
480  Error("InsertRow","column %d out of matrix range",coln);
481  return *this;
482  }
483 
484  if (acoln+nr > fNcols || nr < 0) {
485  Error("InsertRow","row length %d out of range",nr);
486  return *this;
487  }
488  }
489 
490  const Int_t off = arown*fNcols+acoln;
491  Element * const elem = GetMatrixArray()+off;
492  memcpy(elem,v,nr*sizeof(Element));
493 
494  return *this;
495 }
496 
497 ////////////////////////////////////////////////////////////////////////////////
498 /// Store in array v, n matrix elements of row rown starting at column coln
499 
500 template<class Element>
501 void TMatrixTBase<Element>::ExtractRow(Int_t rown,Int_t coln,Element *v,Int_t n) const
502 {
503  const Int_t arown = rown-fRowLwb;
504  const Int_t acoln = coln-fColLwb;
505  const Int_t nr = (n > 0) ? n : fNcols;
506 
507  if (gMatrixCheck) {
508  if (arown >= fNrows || arown < 0) {
509  Error("ExtractRow","row %d out of matrix range",rown);
510  return;
511  }
512 
513  if (acoln >= fNcols || acoln < 0) {
514  Error("ExtractRow","column %d out of matrix range",coln);
515  return;
516  }
517 
518  if (acoln+n >= fNcols || nr < 0) {
519  Error("ExtractRow","row length %d out of range",nr);
520  return;
521  }
522  }
523 
524  const Int_t off = arown*fNcols+acoln;
525  const Element * const elem = GetMatrixArray()+off;
526  memcpy(v,elem,nr*sizeof(Element));
527 }
528 
529 ////////////////////////////////////////////////////////////////////////////////
530 /// Shift the row index by adding row_shift and the column index by adding
531 /// col_shift, respectively. So [rowLwb..rowUpb][colLwb..colUpb] becomes
532 /// [rowLwb+row_shift..rowUpb+row_shift][colLwb+col_shift..colUpb+col_shift]
533 
534 template<class Element>
536 {
537  fRowLwb += row_shift;
538  fColLwb += col_shift;
539 
540  return *this;
541 }
542 
543 ////////////////////////////////////////////////////////////////////////////////
544 /// Set matrix elements to zero
545 
546 template<class Element>
548 {
549  R__ASSERT(IsValid());
550  memset(this->GetMatrixArray(),0,fNelems*sizeof(Element));
551 
552  return *this;
553 }
554 
555 ////////////////////////////////////////////////////////////////////////////////
556 /// Take an absolute value of a matrix, i.e. apply Abs() to each element.
557 
558 template<class Element>
560 {
561  R__ASSERT(IsValid());
562 
563  Element *ep = this->GetMatrixArray();
564  const Element * const fp = ep+fNelems;
565  while (ep < fp) {
566  *ep = TMath::Abs(*ep);
567  ep++;
568  }
569 
570  return *this;
571 }
572 
573 ////////////////////////////////////////////////////////////////////////////////
574 /// Square each element of the matrix.
575 
576 template<class Element>
578 {
579  R__ASSERT(IsValid());
580 
581  Element *ep = this->GetMatrixArray();
582  const Element * const fp = ep+fNelems;
583  while (ep < fp) {
584  *ep = (*ep) * (*ep);
585  ep++;
586  }
587 
588  return *this;
589 }
590 
591 ////////////////////////////////////////////////////////////////////////////////
592 /// Take square root of all elements.
593 
594 template<class Element>
596 {
597  R__ASSERT(IsValid());
598 
599  Element *ep = this->GetMatrixArray();
600  const Element * const fp = ep+fNelems;
601  while (ep < fp) {
602  *ep = TMath::Sqrt(*ep);
603  ep++;
604  }
605 
606  return *this;
607 }
608 
609 ////////////////////////////////////////////////////////////////////////////////
610 /// Make a unit matrix (matrix need not be a square one).
611 
612 template<class Element>
614 {
615  R__ASSERT(IsValid());
616 
617  Element *ep = this->GetMatrixArray();
618  memset(ep,0,fNelems*sizeof(Element));
619  for (Int_t i = fRowLwb; i <= fRowLwb+fNrows-1; i++)
620  for (Int_t j = fColLwb; j <= fColLwb+fNcols-1; j++)
621  *ep++ = (i==j ? 1.0 : 0.0);
622 
623  return *this;
624 }
625 
626 ////////////////////////////////////////////////////////////////////////////////
627 /// option:
628 /// - "D" : b(i,j) = a(i,j)/sqrt(abs*(v(i)*v(j))) (default)
629 /// - else : b(i,j) = a(i,j)*sqrt(abs*(v(i)*v(j))) (default)
630 
631 template<class Element>
633 {
634  R__ASSERT(IsValid());
635  R__ASSERT(v.IsValid());
636 
637  if (gMatrixCheck) {
638  const Int_t nMax = TMath::Max(fNrows,fNcols);
639  if (v.GetNoElements() < nMax) {
640  Error("NormByDiag","vector shorter than matrix diagonal");
641  return *this;
642  }
643  }
644 
645  TString opt(option);
646  opt.ToUpper();
647  const Int_t divide = (opt.Contains("D")) ? 1 : 0;
648 
649  const Element *pV = v.GetMatrixArray();
650  Element *mp = this->GetMatrixArray();
651 
652  if (divide) {
653  for (Int_t irow = 0; irow < fNrows; irow++) {
654  if (pV[irow] != 0.0) {
655  for (Int_t icol = 0; icol < fNcols; icol++) {
656  if (pV[icol] != 0.0) {
657  const Element val = TMath::Sqrt(TMath::Abs(pV[irow]*pV[icol]));
658  *mp++ /= val;
659  } else {
660  Error("NormbyDiag","vector element %d is zero",icol);
661  mp++;
662  }
663  }
664  } else {
665  Error("NormbyDiag","vector element %d is zero",irow);
666  mp += fNcols;
667  }
668  }
669  } else {
670  for (Int_t irow = 0; irow < fNrows; irow++) {
671  for (Int_t icol = 0; icol < fNcols; icol++) {
672  const Element val = TMath::Sqrt(TMath::Abs(pV[irow]*pV[icol]));
673  *mp++ *= val;
674  }
675  }
676  }
677 
678  return *this;
679 }
680 
681 ////////////////////////////////////////////////////////////////////////////////
682 /// Row matrix norm, MAX{ SUM{ |M(i,j)|, over j}, over i}.
683 /// The norm is induced by the infinity vector norm.
684 
685 template<class Element>
687 {
688  R__ASSERT(IsValid());
689 
690  const Element * ep = GetMatrixArray();
691  const Element * const fp = ep+fNelems;
692  Element norm = 0;
693 
694  // Scan the matrix row-after-row
695  while (ep < fp) {
696  Element sum = 0;
697  // Scan a row to compute the sum
698  for (Int_t j = 0; j < fNcols; j++)
699  sum += TMath::Abs(*ep++);
700  norm = TMath::Max(norm,sum);
701  }
702 
703  R__ASSERT(ep == fp);
704 
705  return norm;
706 }
707 
708 ////////////////////////////////////////////////////////////////////////////////
709 /// Column matrix norm, MAX{ SUM{ |M(i,j)|, over i}, over j}.
710 /// The norm is induced by the 1 vector norm.
711 
712 template<class Element>
714 {
715  R__ASSERT(IsValid());
716 
717  const Element * ep = GetMatrixArray();
718  const Element * const fp = ep+fNcols;
719  Element norm = 0;
720 
721  // Scan the matrix col-after-col
722  while (ep < fp) {
723  Element sum = 0;
724  // Scan a col to compute the sum
725  for (Int_t i = 0; i < fNrows; i++,ep += fNcols)
726  sum += TMath::Abs(*ep);
727  ep -= fNelems-1; // Point ep to the beginning of the next col
728  norm = TMath::Max(norm,sum);
729  }
730 
731  R__ASSERT(ep == fp);
732 
733  return norm;
734 }
735 
736 ////////////////////////////////////////////////////////////////////////////////
737 /// Square of the Euclidian norm, SUM{ m(i,j)^2 }.
738 
739 template<class Element>
741 {
742  R__ASSERT(IsValid());
743 
744  const Element * ep = GetMatrixArray();
745  const Element * const fp = ep+fNelems;
746  Element sum = 0;
747 
748  for ( ; ep < fp; ep++)
749  sum += (*ep) * (*ep);
750 
751  return sum;
752 }
753 
754 ////////////////////////////////////////////////////////////////////////////////
755 /// Compute the number of elements != 0.0
756 
757 template<class Element>
759 {
760  R__ASSERT(IsValid());
761 
762  Int_t nr_nonzeros = 0;
763  const Element *ep = this->GetMatrixArray();
764  const Element * const fp = ep+fNelems;
765  while (ep < fp)
766  if (*ep++ != 0.0) nr_nonzeros++;
767 
768  return nr_nonzeros;
769 }
770 
771 ////////////////////////////////////////////////////////////////////////////////
772 /// Compute sum of elements
773 
774 template<class Element>
776 {
777  R__ASSERT(IsValid());
778 
779  Element sum = 0.0;
780  const Element *ep = this->GetMatrixArray();
781  const Element * const fp = ep+fNelems;
782  while (ep < fp)
783  sum += *ep++;
784 
785  return sum;
786 }
787 
788 ////////////////////////////////////////////////////////////////////////////////
789 /// return minimum matrix element value
790 
791 template<class Element>
793 {
794  R__ASSERT(IsValid());
795 
796  const Element * const ep = this->GetMatrixArray();
797  const Int_t index = TMath::LocMin(fNelems,ep);
798  return ep[index];
799 }
800 
801 ////////////////////////////////////////////////////////////////////////////////
802 /// return maximum vector element value
803 
804 template<class Element>
806 {
807  R__ASSERT(IsValid());
808 
809  const Element * const ep = this->GetMatrixArray();
810  const Int_t index = TMath::LocMax(fNelems,ep);
811  return ep[index];
812 }
813 
814 ////////////////////////////////////////////////////////////////////////////////
815 /// Draw this matrix
816 /// The histogram is named "TMatrixT" by default and no title
817 
818 template<class Element>
820 {
821  gROOT->ProcessLine(Form("THistPainter::PaintSpecialObjects((TObject*)0x%lx,\"%s\");",
822  (ULong_t)this, option));
823 }
824 
825 ////////////////////////////////////////////////////////////////////////////////
826 /// Print the matrix as a table of elements.
827 /// By default the format "%11.4g" is used to print one element.
828 /// One can specify an alternative format with eg
829 /// option ="f= %6.2f "
830 
831 template<class Element>
833 {
834  if (!IsValid()) {
835  Error("Print","Matrix is invalid");
836  return;
837  }
838 
839  //build format
840  const char *format = "%11.4g ";
841  if (option) {
842  const char *f = strstr(option,"f=");
843  if (f) format = f+2;
844  }
845  char topbar[100];
846  snprintf(topbar,100,format,123.456789);
847  Int_t nch = strlen(topbar)+1;
848  if (nch > 18) nch = 18;
849  char ftopbar[20];
850  for (Int_t i = 0; i < nch; i++) ftopbar[i] = ' ';
851  Int_t nk = 1 + Int_t(TMath::Log10(fNcols));
852  snprintf(ftopbar+nch/2,20-nch/2,"%s%dd","%",nk);
853  Int_t nch2 = strlen(ftopbar);
854  for (Int_t i = nch2; i < nch; i++) ftopbar[i] = ' ';
855  ftopbar[nch] = '|';
856  ftopbar[nch+1] = 0;
857 
858  printf("\n%dx%d matrix is as follows",fNrows,fNcols);
859 
860  Int_t cols_per_sheet = 5;
861  if (nch <= 8) cols_per_sheet =10;
862  const Int_t ncols = fNcols;
863  const Int_t nrows = fNrows;
864  const Int_t collwb = fColLwb;
865  const Int_t rowlwb = fRowLwb;
866  nk = 5+nch*TMath::Min(cols_per_sheet,fNcols);
867  for (Int_t i = 0; i < nk; i++) topbar[i] = '-';
868  topbar[nk] = 0;
869  for (Int_t sheet_counter = 1; sheet_counter <= ncols; sheet_counter += cols_per_sheet) {
870  printf("\n\n |");
871  for (Int_t j = sheet_counter; j < sheet_counter+cols_per_sheet && j <= ncols; j++)
872  printf(ftopbar,j+collwb-1);
873  printf("\n%s\n",topbar);
874  if (fNelems <= 0) continue;
875  for (Int_t i = 1; i <= nrows; i++) {
876  printf("%4d |",i+rowlwb-1);
877  for (Int_t j = sheet_counter; j < sheet_counter+cols_per_sheet && j <= ncols; j++)
878  printf(format,(*this)(i+rowlwb-1,j+collwb-1));
879  printf("\n");
880  }
881  }
882  printf("\n");
883 }
884 
885 ////////////////////////////////////////////////////////////////////////////////
886 /// Are all matrix elements equal to val?
887 
888 template<class Element>
890 {
891  R__ASSERT(IsValid());
892 
893  if (val == 0. && fNelems == 0)
894  return kTRUE;
895 
896  const Element * ep = GetMatrixArray();
897  const Element * const fp = ep+fNelems;
898  for (; ep < fp; ep++)
899  if (!(*ep == val))
900  return kFALSE;
901 
902  return kTRUE;
903 }
904 
905 ////////////////////////////////////////////////////////////////////////////////
906 /// Are all matrix elements not equal to val?
907 
908 template<class Element>
910 {
911  R__ASSERT(IsValid());
912 
913  if (val == 0. && fNelems == 0)
914  return kFALSE;
915 
916  const Element * ep = GetMatrixArray();
917  const Element * const fp = ep+fNelems;
918  for (; ep < fp; ep++)
919  if (!(*ep != val))
920  return kFALSE;
921 
922  return kTRUE;
923 }
924 
925 ////////////////////////////////////////////////////////////////////////////////
926 /// Are all matrix elements < val?
927 
928 template<class Element>
930 {
931  R__ASSERT(IsValid());
932 
933  const Element * ep = GetMatrixArray();
934  const Element * const fp = ep+fNelems;
935  for (; ep < fp; ep++)
936  if (!(*ep < val))
937  return kFALSE;
938 
939  return kTRUE;
940 }
941 
942 ////////////////////////////////////////////////////////////////////////////////
943 /// Are all matrix elements <= val?
944 
945 template<class Element>
947 {
948  R__ASSERT(IsValid());
949 
950  const Element * ep = GetMatrixArray();
951  const Element * const fp = ep+fNelems;
952  for (; ep < fp; ep++)
953  if (!(*ep <= val))
954  return kFALSE;
955 
956  return kTRUE;
957 }
958 
959 ////////////////////////////////////////////////////////////////////////////////
960 /// Are all matrix elements > val?
961 
962 template<class Element>
964 {
965  R__ASSERT(IsValid());
966 
967  const Element * ep = GetMatrixArray();
968  const Element * const fp = ep+fNelems;
969  for (; ep < fp; ep++)
970  if (!(*ep > val))
971  return kFALSE;
972 
973  return kTRUE;
974 }
975 
976 ////////////////////////////////////////////////////////////////////////////////
977 /// Are all matrix elements >= val?
978 
979 template<class Element>
981 {
982  R__ASSERT(IsValid());
983 
984  const Element * ep = GetMatrixArray();
985  const Element * const fp = ep+fNelems;
986  for (; ep < fp; ep++)
987  if (!(*ep >= val))
988  return kFALSE;
989 
990  return kTRUE;
991 }
992 
993 ////////////////////////////////////////////////////////////////////////////////
994 /// Apply action to each matrix element
995 
996 template<class Element>
998 {
999  R__ASSERT(IsValid());
1000 
1001  Element *ep = this->GetMatrixArray();
1002  const Element * const ep_last = ep+fNelems;
1003  while (ep < ep_last)
1004  action.Operation(*ep++);
1005 
1006  return *this;
1007 }
1008 
1009 ////////////////////////////////////////////////////////////////////////////////
1010 /// Apply action to each element of the matrix. To action the location
1011 /// of the current element is passed.
1012 
1013 template<class Element>
1015 {
1016  R__ASSERT(IsValid());
1017 
1018  Element *ep = this->GetMatrixArray();
1019  for (action.fI = fRowLwb; action.fI < fRowLwb+fNrows; action.fI++)
1020  for (action.fJ = fColLwb; action.fJ < fColLwb+fNcols; action.fJ++)
1021  action.Operation(*ep++);
1022 
1023  R__ASSERT(ep == this->GetMatrixArray()+fNelems);
1024 
1025  return *this;
1026 }
1027 
1028 ////////////////////////////////////////////////////////////////////////////////
1029 /// Randomize matrix element values
1030 
1031 template<class Element>
1033 {
1034  R__ASSERT(IsValid());
1035 
1036  const Element scale = beta-alpha;
1037  const Element shift = alpha/scale;
1038 
1039  Element * ep = GetMatrixArray();
1040  const Element * const fp = ep+fNelems;
1041  while (ep < fp)
1042  *ep++ = scale*(Drand(seed)+shift);
1043 
1044  return *this;
1045 }
1046 
1047 ////////////////////////////////////////////////////////////////////////////////
1048 /// Check to see if two matrices are identical.
1049 
1050 template<class Element>
1052 {
1053  if (!AreCompatible(m1,m2)) return kFALSE;
1054  return (memcmp(m1.GetMatrixArray(),m2.GetMatrixArray(),
1055  m1.GetNoElements()*sizeof(Element)) == 0);
1056 }
1057 
1058 ////////////////////////////////////////////////////////////////////////////////
1059 /// Square of the Euclidian norm of the difference between two matrices.
1060 
1061 template<class Element>
1063 {
1064  if (gMatrixCheck && !AreCompatible(m1,m2)) {
1065  ::Error("E2Norm","matrices not compatible");
1066  return -1.0;
1067  }
1068 
1069  const Element * mp1 = m1.GetMatrixArray();
1070  const Element * mp2 = m2.GetMatrixArray();
1071  const Element * const fmp1 = mp1+m1.GetNoElements();
1072 
1073  Element sum = 0.0;
1074  for (; mp1 < fmp1; mp1++, mp2++)
1075  sum += (*mp1 - *mp2)*(*mp1 - *mp2);
1076 
1077  return sum;
1078 }
1079 
1080 ////////////////////////////////////////////////////////////////////////////////
1081 /// Check that matrice sm1 and m2 areboth valid and have identical shapes .
1082 
1083 template<class Element1,class Element2>
1085 {
1086  if (!m1.IsValid()) {
1087  if (verbose)
1088  ::Error("AreCompatible", "matrix 1 not valid");
1089  return kFALSE;
1090  }
1091  if (!m2.IsValid()) {
1092  if (verbose)
1093  ::Error("AreCompatible", "matrix 2 not valid");
1094  return kFALSE;
1095  }
1096 
1097  if (m1.GetNrows() != m2.GetNrows() || m1.GetNcols() != m2.GetNcols() ||
1098  m1.GetRowLwb() != m2.GetRowLwb() || m1.GetColLwb() != m2.GetColLwb()) {
1099  if (verbose)
1100  ::Error("AreCompatible", "matrices 1 and 2 not compatible");
1101  return kFALSE;
1102  }
1103 
1104  return kTRUE;
1105 }
1106 
1107 ////////////////////////////////////////////////////////////////////////////////
1108 /// Compare two matrices and print out the result of the comparison.
1109 
1110 template<class Element>
1112 {
1113  if (!AreCompatible(m1,m2)) {
1114  Error("Compare(const TMatrixTBase<Element> &,const TMatrixTBase<Element> &)","matrices are incompatible");
1115  return;
1116  }
1117 
1118  printf("\n\nComparison of two TMatrices:\n");
1119 
1120  Element norm1 = 0; // Norm of the Matrices
1121  Element norm2 = 0;
1122  Element ndiff = 0; // Norm of the difference
1123  Int_t imax = 0; // For the elements that differ most
1124  Int_t jmax = 0;
1125  Element difmax = -1;
1126 
1127  for (Int_t i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++) {
1128  for (Int_t j = m1.GetColLwb(); j < m1.GetColUpb(); j++) {
1129  const Element mv1 = m1(i,j);
1130  const Element mv2 = m2(i,j);
1131  const Element diff = TMath::Abs(mv1-mv2);
1132 
1133  if (diff > difmax) {
1134  difmax = diff;
1135  imax = i;
1136  jmax = j;
1137  }
1138  norm1 += TMath::Abs(mv1);
1139  norm2 += TMath::Abs(mv2);
1140  ndiff += TMath::Abs(diff);
1141  }
1142  }
1143 
1144  printf("\nMaximal discrepancy \t\t%g", difmax);
1145  printf("\n occured at the point\t\t(%d,%d)",imax,jmax);
1146  const Element mv1 = m1(imax,jmax);
1147  const Element mv2 = m2(imax,jmax);
1148  printf("\n Matrix 1 element is \t\t%g", mv1);
1149  printf("\n Matrix 2 element is \t\t%g", mv2);
1150  printf("\n Absolute error v2[i]-v1[i]\t\t%g", mv2-mv1);
1151  printf("\n Relative error\t\t\t\t%g\n",
1152  (mv2-mv1)/TMath::Max(TMath::Abs(mv2+mv1)/2,(Element)1e-7));
1153 
1154  printf("\n||Matrix 1|| \t\t\t%g", norm1);
1155  printf("\n||Matrix 2|| \t\t\t%g", norm2);
1156  printf("\n||Matrix1-Matrix2||\t\t\t\t%g", ndiff);
1157  printf("\n||Matrix1-Matrix2||/sqrt(||Matrix1|| ||Matrix2||)\t%g\n\n",
1158  ndiff/TMath::Max(TMath::Sqrt(norm1*norm2),1e-7));
1159 }
1160 
1161 ////////////////////////////////////////////////////////////////////////////////
1162 /// Validate that all elements of matrix have value val within maxDevAllow.
1163 
1164 template<class Element>
1165 Bool_t VerifyMatrixValue(const TMatrixTBase<Element> &m,Element val,Int_t verbose,Element maxDevAllow)
1166 {
1167  R__ASSERT(m.IsValid());
1168 
1169  if (m == 0)
1170  return kTRUE;
1171 
1172  Int_t imax = 0;
1173  Int_t jmax = 0;
1174  Element maxDevObs = 0;
1175 
1176  if (TMath::Abs(maxDevAllow) <= 0.0)
1177  maxDevAllow = std::numeric_limits<Element>::epsilon();
1178 
1179  for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
1180  for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++) {
1181  const Element dev = TMath::Abs(m(i,j)-val);
1182  if (dev > maxDevObs) {
1183  imax = i;
1184  jmax = j;
1185  maxDevObs = dev;
1186  }
1187  }
1188  }
1189 
1190  if (maxDevObs == 0)
1191  return kTRUE;
1192 
1193  if (verbose) {
1194  printf("Largest dev for (%d,%d); dev = |%g - %g| = %g\n",imax,jmax,m(imax,jmax),val,maxDevObs);
1195  if(maxDevObs > maxDevAllow)
1196  Error("VerifyElementValue","Deviation > %g\n",maxDevAllow);
1197  }
1198 
1199  if(maxDevObs > maxDevAllow)
1200  return kFALSE;
1201  return kTRUE;
1202 }
1203 
1204 ////////////////////////////////////////////////////////////////////////////////
1205 /// Verify that elements of the two matrices are equal within MaxDevAllow .
1206 
1207 template<class Element>
1209  Element maxDevAllow)
1210 {
1211  if (!AreCompatible(m1,m2,verbose))
1212  return kFALSE;
1213 
1214  if (m1 == 0 && m2 == 0)
1215  return kTRUE;
1216 
1217  Int_t imax = 0;
1218  Int_t jmax = 0;
1219  Element maxDevObs = 0;
1220 
1221  if (TMath::Abs(maxDevAllow) <= 0.0)
1222  maxDevAllow = std::numeric_limits<Element>::epsilon();
1223 
1224  for (Int_t i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++) {
1225  for (Int_t j = m1.GetColLwb(); j <= m1.GetColUpb(); j++) {
1226  const Element dev = TMath::Abs(m1(i,j)-m2(i,j));
1227  if (dev > maxDevObs) {
1228  imax = i;
1229  jmax = j;
1230  maxDevObs = dev;
1231  }
1232  }
1233  }
1234 
1235  if (maxDevObs == 0)
1236  return kTRUE;
1237 
1238  if (verbose) {
1239  printf("Largest dev for (%d,%d); dev = |%g - %g| = %g\n",
1240  imax,jmax,m1(imax,jmax),m2(imax,jmax),maxDevObs);
1241  if (maxDevObs > maxDevAllow)
1242  Error("VerifyMatrixValue","Deviation > %g\n",maxDevAllow);
1243  }
1244 
1245  if (maxDevObs > maxDevAllow)
1246  return kFALSE;
1247  return kTRUE;
1248 }
1249 
1250 ////////////////////////////////////////////////////////////////////////////////
1251 /// Stream an object of class TMatrixTBase<Element>.
1252 
1253 template<class Element>
1255 {
1256  if (R__b.IsReading()) {
1257  UInt_t R__s, R__c;
1258  Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
1259  if (R__v > 1) {
1260  R__b.ReadClassBuffer(TMatrixTBase<Element>::Class(),this,R__v,R__s,R__c);
1261  } else {
1262  Error("TMatrixTBase<Element>::Streamer","Unknown version number: %d",R__v);
1263  R__ASSERT(0);
1264  }
1265  if (R__v < 4) MakeValid();
1266  } else {
1268  }
1269 }
1270 
1271 // trick to return a reference to nan in operator(i,j_ when i,j are outside of range
1272 template<class Element>
1273 struct nan_value_t {
1274  static Element gNanValue;
1275 };
1276 template<>
1277 Double_t nan_value_t<Double_t>::gNanValue = std::numeric_limits<Double_t>::quiet_NaN();
1278 template<>
1279 Float_t nan_value_t<Float_t>::gNanValue = std::numeric_limits<Float_t>::quiet_NaN();
1280 
1281 template<class Element>
1283 {
1284  return nan_value_t<Element>::gNanValue;
1285 }
1286 
1287 
1288 template class TMatrixTBase<Float_t>;
1289 
1290 template Bool_t operator== <Float_t>(const TMatrixFBase &m1,const TMatrixFBase &m2);
1291 template Float_t E2Norm <Float_t>(const TMatrixFBase &m1,const TMatrixFBase &m2);
1293  (const TMatrixFBase &m1,const TMatrixFBase &m2,Int_t verbose);
1295  (const TMatrixFBase &m1,const TMatrixDBase &m2,Int_t verbose);
1296 template void Compare <Float_t>(const TMatrixFBase &m1,const TMatrixFBase &m2);
1297 template Bool_t VerifyMatrixValue <Float_t>(const TMatrixFBase &m,Float_t val,Int_t verbose,Float_t maxDevAllow);
1300  Int_t verbose,Float_t maxDevAllowN);
1302 
1303 template class TMatrixTBase<Double_t>;
1304 
1305 template Bool_t operator== <Double_t>(const TMatrixDBase &m1,const TMatrixDBase &m2);
1306 template Double_t E2Norm <Double_t>(const TMatrixDBase &m1,const TMatrixDBase &m2);
1308  (const TMatrixDBase &m1,const TMatrixDBase &m2,Int_t verbose);
1310  (const TMatrixDBase &m1,const TMatrixFBase &m2,Int_t verbose);
1311 template void Compare <Double_t>(const TMatrixDBase &m1,const TMatrixDBase &m2);
1315  Int_t verbose,Double_t maxDevAllow);
Bool_t VerifyMatrixIdentity(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2, Int_t verbose, Element maxDevAllow)
Verify that elements of the two matrices are equal within MaxDevAllow .
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
virtual const Element * GetMatrixArray() const =0
virtual TMatrixTBase< Element > & Sqrt()
Take square root of all elements.
virtual Bool_t IsSymmetric() const
Check whether matrix is symmetric.
Bool_t IsReading() const
Definition: TBuffer.h:85
virtual TMatrixTBase< Element > & Shift(Int_t row_shift, Int_t col_shift)
Shift the row index by adding row_shift and the column index by adding col_shift, respectively...
virtual Int_t NonZeros() const
Compute the number of elements != 0.0.
static long int sum(long int i)
Definition: Factory.cxx:2276
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
template Bool_t AreCompatible< Float_t, Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2, Int_t verbose)
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:122
auto * m
Definition: textangle.C:8
Long64_t LocMax(Long64_t n, const T *a)
Return index of array with the maximum element.
Definition: TMath.h:990
virtual TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *option="")
Copy array data to matrix .
short Version_t
Definition: RtypesCore.h:61
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
Bool_t operator!=(Element val) const
Are all matrix elements not equal to val?
float Float_t
Definition: RtypesCore.h:53
const char Option_t
Definition: RtypesCore.h:62
Int_t GetNcols() const
Definition: TMatrixTBase.h:127
TVectorT.
Definition: TMatrixTBase.h:79
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1138
Buffer base class used for serializing objects.
Definition: TBuffer.h:42
#define R__ASSERT(e)
Definition: TError.h:96
#define gROOT
Definition: TROOT.h:415
Basic string class.
Definition: TString.h:131
#define f(i)
Definition: RSha256.hxx:104
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
virtual TMatrixTBase< Element > & Sqr()
Square each element of the matrix.
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
Bool_t operator>(Element val) const
Are all matrix elements > val?
void Compare(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Compare two matrices and print out the result of the comparison.
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
double beta(double x, double y)
Calculates the beta function.
template Double_t E2Norm< Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2)
template Bool_t VerifyMatrixValue< Float_t >(const TMatrixFBase &m, Float_t val, Int_t verbose, Float_t maxDevAllow)
Int_t GetNoElements() const
Definition: TMatrixTBase.h:128
Bool_t operator<(Element val) const
Are all matrix elements < val?
static void DoubleLexSort(Int_t n, Int_t *first, Int_t *second, Element *data)
default kTRUE, when Use array kFALSE
#define templateClassImp(name)
Definition: Rtypes.h:409
Int_t GetColLwb() const
Definition: TMatrixTBase.h:125
virtual void Operation(Element &element) const =0
Element E2Norm(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Square of the Euclidian norm of the difference between two matrices.
Double_t Log10(Double_t x)
Definition: TMath.h:754
Bool_t operator>=(Element val) const
Are all matrix elements >= val?
template Bool_t AreCompatible< Double_t, Float_t >(const TMatrixDBase &m1, const TMatrixFBase &m2, Int_t verbose)
void Error(const char *location, const char *msgfmt,...)
Bool_t operator==(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Check to see if two matrices are identical.
template Bool_t VerifyMatrixIdentity< Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2, Int_t verbose, Double_t maxDevAllow)
virtual TMatrixTBase< Element > & Apply(const TElementActionT< Element > &action)
Apply action to each matrix element.
virtual TMatrixTBase< Element > & NormByDiag(const TVectorT< Element > &v, Option_t *option="D")
option:
static constexpr double m2
template Bool_t AreCompatible< Double_t, Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2, Int_t verbose)
virtual void ExtractRow(Int_t row, Int_t col, Element *v, Int_t n=-1) const
Store in array v, n matrix elements of row rown starting at column coln.
virtual TMatrixTBase< Element > & InsertRow(Int_t row, Int_t col, const Element *v, Int_t n=-1)
Copy n elements from array v to row rown starting at column coln.
unsigned int UInt_t
Definition: RtypesCore.h:42
char * Form(const char *fmt,...)
virtual void GetMatrix2Array(Element *data, Option_t *option="") const
Copy matrix data to array .
virtual Element Min() const
return minimum matrix element value
virtual TMatrixTBase< Element > & Randomize(Element alpha, Element beta, Double_t &seed)
Randomize matrix element values.
template Bool_t AreCompatible< Float_t, Double_t >(const TMatrixFBase &m1, const TMatrixDBase &m2, Int_t verbose)
Bool_t operator<=(Element val) const
Are all matrix elements <= val?
Int_t GetColUpb() const
Definition: TMatrixTBase.h:126
Linear Algebra Package.
REAL epsilon
Definition: triangle.c:617
Int_t GetNrows() const
Definition: TMatrixTBase.h:124
const Bool_t kFALSE
Definition: RtypesCore.h:88
virtual Element RowNorm() const
Row matrix norm, MAX{ SUM{ |M(i,j)|, over j}, over i}.
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
virtual Element E2Norm() const
Square of the Euclidian norm, SUM{ m(i,j)^2 }.
template Bool_t VerifyMatrixIdentity< Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2, Int_t verbose, Float_t maxDevAllowN)
virtual Element Sum() const
Compute sum of elements.
template void Compare< Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2)
double Double_t
Definition: RtypesCore.h:55
template void Compare< Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2)
unsigned long ULong_t
Definition: RtypesCore.h:51
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Bool_t AreCompatible(const TMatrixTBase< Element1 > &m1, const TMatrixTBase< Element2 > &m2, Int_t verbose)
Check that matrice sm1 and m2 areboth valid and have identical shapes .
virtual void Operation(Element &element) const =0
void Print(Option_t *name="") const
Print the matrix as a table of elements.
void Draw(Option_t *option="")
Draw this matrix The histogram is named "TMatrixT" by default and no title.
Int_t gMatrixCheck
TMatrixTBase.
static void IndexedLexSort(Int_t n, Int_t *first, Int_t swapFirst, Int_t *second, Int_t swapSecond, Int_t *index)
Lexical sort on array data using indices first and second.
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
Bool_t operator==(Element val) const
Are all matrix elements equal to val?
#define snprintf
Definition: civetweb.c:1540
static Element & NaNValue()
Int_t GetRowUpb() const
Definition: TMatrixTBase.h:123
Definition: first.py:1
virtual Element ColNorm() const
Column matrix norm, MAX{ SUM{ |M(i,j)|, over i}, over j}.
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
Long64_t LocMin(Long64_t n, const T *a)
Return index of array with the minimum element.
Definition: TMath.h:962
const Bool_t kTRUE
Definition: RtypesCore.h:87
template Bool_t VerifyMatrixValue< Double_t >(const TMatrixDBase &m, Double_t val, Int_t verbose, Double_t maxDevAllow)
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
virtual Element Max() const
return maximum vector element value
Bool_t VerifyMatrixValue(const TMatrixTBase< Element > &m, Element val, Int_t verbose, Element maxDevAllow)
Validate that all elements of matrix have value val within maxDevAllow.
template Float_t E2Norm< Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2)
const Int_t n
Definition: legend1.C:16
Bool_t IsValid() const
Definition: TMatrixTBase.h:147
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
virtual TMatrixTBase< Element > & Abs()
Take an absolute value of a matrix, i.e. apply Abs() to each element.
static constexpr double second