ROOT   6.08/07 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 /*************************************************************************
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
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>
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 {
1257  UInt_t R__s, R__c;
1258  Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
1259  if (R__v > 1) {
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);
1299 template Bool_t VerifyMatrixIdentity<Float_t>(const TMatrixFBase &m1,const TMatrixFBase &m2,
1300  Int_t verbose,Float_t maxDevAllowN);
1301 template Bool_t VerifyMatrixIdentity<Float_t>(const TMatrixFBase &m1,const TMatrixFBase &m2);
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);
1314 template Bool_t VerifyMatrixIdentity<Double_t>(const TMatrixDBase &m1,const TMatrixDBase &m2,
1315  Int_t verbose,Double_t maxDevAllow);
1316 template Bool_t VerifyMatrixIdentity<Double_t>(const TMatrixDBase &m1,const TMatrixDBase &m2);
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.
Definition: TBuffer.h:83
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:1786
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:132
Long64_t LocMax(Long64_t n, const T *a)
Definition: TMath.h:711
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:137
TVectorT.
Definition: TMatrixTBase.h:89
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1102
Buffer base class used for serializing objects.
Definition: TBuffer.h:42
#define R__ASSERT(e)
Definition: TError.h:98
#define gROOT
Definition: TROOT.h:364
Basic string class.
Definition: TString.h:137
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:170
virtual TMatrixTBase< Element > & Sqr()
Square each element of the matrix.
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
const Bool_t kFALSE
Definition: Rtypes.h:92
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.
static std::string format(double x, double y, int digits, int width)
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
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:138
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:325
Int_t GetColLwb() const
Definition: TMatrixTBase.h:135
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:529
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)
Element * GetMatrixArray()
Definition: TVectorT.h:84
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:
SVector< double, 2 > v
Definition: Dict.h:5
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
Bool_t IsValid() const
Definition: TVectorT.h:89
TMarker * m
Definition: textangle.C:8
bool verbose
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.
Int_t GetNoElements() const
Definition: TVectorT.h:82
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:136
Linear Algebra Package.
REAL epsilon
Definition: triangle.c:617
Int_t GetNrows() const
Definition: TMatrixTBase.h:134
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 }.
double f(double x)
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:567
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:202
Bool_t operator==(Element val) const
Are all matrix elements equal to val?
#define snprintf
Definition: civetweb.c:822
static Element & NaNValue()
Int_t GetRowUpb() const
Definition: TMatrixTBase.h:133
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:464
Long64_t LocMin(Long64_t n, const T *a)
Definition: TMath.h:682
const Bool_t kTRUE
Definition: Rtypes.h:91
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40
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:157
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.