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