Logo ROOT   6.18/05
Reference Guide
TMatrixTBase.cxx
Go to the documentation of this file.
1// @(#)root/matrix:$Id: 2d00df45ce4c38c7ea0930d6b520cbf4cfb9152e $
2// Authors: Fons Rademakers, Eddy Offermann Nov 2003
3
4/*************************************************************************
5 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/** \class TMatrixTBase
13 \ingroup Matrix
14 \brief Linear Algebra Package
15
16 ### Linear Algebra Package
17
18 The present package implements all the basic algorithms dealing
19 with vectors, matrices, matrix columns, rows, diagonals, etc.
20 In addition eigen-Vector analysis and several matrix decomposition
21 have been added (LU,QRH,Cholesky,Bunch-Kaufman and SVD) .
22 The decompositions are used in matrix inversion, equation solving.
23
24 For a dense matrix, elements are arranged in memory in a ROW-wise
25 fashion . For (n x m) matrices where n*m <=kSizeMax (=25 currently)
26 storage space is available on the stack, thus avoiding expensive
27 allocation/deallocation of heap space . However, this introduces of
28 course kSizeMax overhead for each matrix object . If this is an
29 issue recompile with a new appropriate value (>=0) for kSizeMax
30
31 Sparse matrices are also stored in row-wise fashion but additional
32 row/column information is stored, see TMatrixTSparse source for
33 additional details .
34
35 Another way to assign and store matrix data is through Use
36 see for instance stressLinear.cxx file .
37
38 Unless otherwise specified, matrix and vector indices always start
39 with 0, spanning up to the specified limit-1. However, there are
40 constructors to which one can specify aribtrary lower and upper
41 bounds, e.g. TMatrixD m(1,10,1,5) defines a matrix that ranges
42 from 1..10, 1..5 (a(1,1)..a(10,5)).
43
44 The present package provides all facilities to completely AVOID
45 returning matrices. Use "TMatrixD A(TMatrixD::kTransposed,B);"
46 and other fancy constructors as much as possible. If one really needs
47 to return a matrix, return a TMatrixTLazy object instead. The
48 conversion is completely transparent to the end user, e.g.
49 "TMatrixT m = THaarMatrixT(5);" and _is_ efficient.
50
51 Since TMatrixT et al. are fully integrated in ROOT, they of course
52 can be stored in a ROOT database.
53
54 For usage examples see $ROOTSYS/test/stressLinear.cxx
55
56 ### Acknowledgements
57
58 1. Oleg E. Kiselyov
59 First implementations were based on the his code . We have diverged
60 quite a bit since then but the ideas/code for lazy matrix and
61 "nested function" are 100% his .
62 You can see him and his code in action at http://okmij.org/ftp
63 2. Chris R. Birchenhall,
64 We adapted his idea of the implementation for the decomposition
65 classes instead of our messy installation of matrix inversion
66 His installation of matrix condition number, using an iterative
67 scheme using the Hage algorithm is worth looking at !
68 Chris has a nice writeup (matdoc.ps) on his matrix classes at
69 ftp://ftp.mcc.ac.uk/pub/matclass/
70 3. Mark Fischler and Steven Haywood of CLHEP
71 They did the slave labor of spelling out all sub-determinants
72 for Cramer inversion of (4x4),(5x5) and (6x6) matrices
73 The stack storage for small matrices was also taken from them
74 4. Roldan Pozo of TNT (http://math.nist.gov/tnt/)
75 He converted the EISPACK routines for the eigen-vector analysis to
76 C++ . We started with his implementation
77 5. Siegmund Brandt (http://siux00.physik.uni-siegen.de/~brandt/datan
78 We adapted his (very-well) documented SVD routines
79
80### How to efficiently use this package
81
82#### 1. Never return complex objects (matrices or vectors)
83 Danger: For example, when the following snippet:
84~~~
85 TMatrixD foo(int n)
86 {
87 TMatrixD foom(n,n); fill_in(foom); return foom;
88 }
89 TMatrixD m = foo(5);
90~~~
91 runs, it constructs matrix foo:foom, copies it onto stack as a
92 return value and destroys foo:foom. Return value (a matrix)
93 from foo() is then copied over to m (via a copy constructor),
94 and the return value is destroyed. So, the matrix constructor is
95 called 3 times and the destructor 2 times. For big matrices,
96 the cost of multiple constructing/copying/destroying of objects
97 may be very large. *Some* optimized compilers can cut down on 1
98 copying/destroying, but still it leaves at least two calls to
99 the constructor. Note, TMatrixDLazy (see below) can construct
100 TMatrixD m "inplace", with only a _single_ call to the
101 constructor.
102
103#### 2. Use "two-address instructions"
104~~~
105 "void TMatrixD::operator += (const TMatrixD &B);"
106~~~
107 as much as possible.
108 That is, to add two matrices, it's much more efficient to write
109~~~
110 A += B;
111~~~
112 than
113~~~
114 TMatrixD C = A + B;
115~~~
116 (if both operand should be preserved,
117 TMatrixD C = A; C += B;
118 is still better).
119
120#### 3. Use glorified constructors when returning of an object seems inevitable:
121~~~
122 "TMatrixD A(TMatrixD::kTransposed,B);"
123 "TMatrixD C(A,TMatrixD::kTransposeMult,B);"
124~~~
125
126 like in the following snippet (from $ROOTSYS/test/vmatrix.cxx)
127 that verifies that for an orthogonal matrix T, T'T = TT' = E.
128
129~~~
130 TMatrixD haar = THaarMatrixD(5);
131 TMatrixD unit(TMatrixD::kUnit,haar);
132 TMatrixD haar_t(TMatrixD::kTransposed,haar);
133 TMatrixD hth(haar,TMatrixD::kTransposeMult,haar);
134 TMatrixD hht(haar,TMatrixD::kMult,haar_t);
135 TMatrixD hht1 = haar; hht1 *= haar_t;
136 VerifyMatrixIdentity(unit,hth);
137 VerifyMatrixIdentity(unit,hht);
138 VerifyMatrixIdentity(unit,hht1);
139~~~
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)
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~~~
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.
175 Another example is applying of a simple function to each matrix element:
177~~~
178 void foo(TMatrixD &m,TMatrixD &m1)
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);
191~~~
193 Validation code $ROOTSYS/test/vmatrix.cxx and vvector.cxx contain
194 a few more examples of that kind.
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
237template<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
276template<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
383template<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
411template<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
440template<class Element>
441void TMatrixTBase<Element>::GetMatrix2Array(Element *data,Option_t *option) const
442{
443 R__ASSERT(IsValid());
444
445 TString opt = option;
446 opt.ToUpper();
447
448 const Element * const elem = GetMatrixArray();
449 if (opt.Contains("F")) {
450 for (Int_t irow = 0; irow < fNrows; irow++) {
451 const Int_t off1 = irow*fNcols;
452 Int_t off2 = 0;
453 for (Int_t icol = 0; icol < fNcols; icol++) {
454 data[off2+irow] = elem[off1+icol];
455 off2 += fNrows;
456 }
457 }
458 }
459 else
460 memcpy(data,elem,fNelems*sizeof(Element));
461}
462
463////////////////////////////////////////////////////////////////////////////////
464/// Copy n elements from array v to row rown starting at column coln
465
466template<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
500template<class Element>
501void 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
534template<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
546template<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
558template<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
576template<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
594template<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
612template<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
631template<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
685template<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
712template<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
739template<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
757template<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
774template<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
791template<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
804template<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
818template<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
831template<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
888template<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
908template<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
928template<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
945template<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
962template<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
979template<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
996template<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
1013template<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
1031template<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
1050template<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
1061template<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
1083template<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
1110template<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
1164template<class Element>
1165Bool_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)
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
1207template<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)
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
1253template<class Element>
1255{
1256 if (R__b.IsReading()) {
1257 UInt_t R__s, R__c;
1258 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
1259 if (R__v > 1) {
1260 R__b.ReadClassBuffer(TMatrixTBase<Element>::Class(),this,R__v,R__s,R__c);
1261 } else {
1262 Error("TMatrixTBase<Element>::Streamer","Unknown version number: %d",R__v);
1263 R__ASSERT(0);
1264 }
1265 if (R__v < 4) MakeValid();
1266 } else {
1268 }
1269}
1270
1271// trick to return a reference to nan in operator(i,j_ when i,j are outside of range
1272template<class Element>
1273struct nan_value_t {
1274 static Element gNanValue;
1275};
1276template<>
1277Double_t nan_value_t<Double_t>::gNanValue = std::numeric_limits<Double_t>::quiet_NaN();
1278template<>
1279Float_t nan_value_t<Float_t>::gNanValue = std::numeric_limits<Float_t>::quiet_NaN();
1280
1281template<class Element>
1283{
1284 return nan_value_t<Element>::gNanValue;
1285}
1286
1287
1288template class TMatrixTBase<Float_t>;
1289
1290template Bool_t operator== <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);
1296template void Compare <Float_t>(const TMatrixFBase &m1,const TMatrixFBase &m2);
1300 Int_t verbose,Float_t maxDevAllowN);
1302
1303template class TMatrixTBase<Double_t>;
1304
1305template Bool_t operator== <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);
1311template void Compare <Double_t>(const TMatrixDBase &m1,const TMatrixDBase &m2);
1315 Int_t verbose,Double_t maxDevAllow);
SVector< double, 2 > v
Definition: Dict.h:5
#define f(i)
Definition: RSha256.hxx:104
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:41
short Version_t
Definition: RtypesCore.h:61
unsigned int UInt_t
Definition: RtypesCore.h:42
const Bool_t kFALSE
Definition: RtypesCore.h:88
unsigned long ULong_t
Definition: RtypesCore.h:51
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
float Float_t
Definition: RtypesCore.h:53
const Bool_t kTRUE
Definition: RtypesCore.h:87
const char Option_t
Definition: RtypesCore.h:62
#define templateClassImp(name)
Definition: Rtypes.h:409
#define R__ASSERT(e)
Definition: TError.h:96
void Error(const char *location, const char *msgfmt,...)
template void Compare< Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2)
template Bool_t VerifyMatrixValue< Double_t >(const TMatrixDBase &m, Double_t val, Int_t verbose, Double_t maxDevAllow)
template void Compare< 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 .
template Bool_t AreCompatible< Float_t, Double_t >(const TMatrixFBase &m1, const TMatrixDBase &m2, Int_t verbose)
template Bool_t AreCompatible< Float_t, Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2, Int_t verbose)
void Compare(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Compare two matrices and print out the result of the comparison.
Element E2Norm(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Square of the Euclidian norm of the difference between two matrices.
template Bool_t VerifyMatrixIdentity< Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2, Int_t verbose, Double_t maxDevAllow)
Int_t gMatrixCheck
TMatrixTBase.
template Double_t E2Norm< Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2)
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)
template Bool_t AreCompatible< Double_t, Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2, Int_t verbose)
template Bool_t VerifyMatrixIdentity< Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2, Int_t verbose, Float_t maxDevAllowN)
template Bool_t AreCompatible< Double_t, Float_t >(const TMatrixDBase &m1, const TMatrixFBase &m2, Int_t verbose)
Bool_t operator==(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Check to see if two matrices are identical.
template Bool_t VerifyMatrixValue< Float_t >(const TMatrixFBase &m, Float_t val, Int_t verbose, Float_t maxDevAllow)
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 .
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
#define gROOT
Definition: TROOT.h:414
char * Form(const char *fmt,...)
#define snprintf
Definition: civetweb.c:1540
Buffer base class used for serializing objects.
Definition: TBuffer.h:42
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
Bool_t IsReading() const
Definition: TBuffer.h:85
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
virtual void Operation(Element &element) const =0
virtual void Operation(Element &element) const =0
Linear Algebra Package.
Definition: TMatrixTBase.h:85
virtual Element Sum() const
Compute sum of elements.
virtual Element RowNorm() const
Row matrix norm, MAX{ SUM{ |M(i,j)|, over j}, over i}.
virtual Element ColNorm() const
Column matrix norm, MAX{ SUM{ |M(i,j)|, over i}, over j}.
void Draw(Option_t *option="")
Draw this matrix The histogram is named "TMatrixT" by default and no title.
virtual TMatrixTBase< Element > & Randomize(Element alpha, Element beta, Double_t &seed)
Randomize matrix element values.
virtual TMatrixTBase< Element > & Sqr()
Square each element of the matrix.
virtual const Element * GetMatrixArray() const =0
Int_t GetNrows() const
Definition: TMatrixTBase.h:124
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
void Print(Option_t *name="") const
Print the matrix as a table of elements.
Bool_t operator!=(Element val) const
Are all matrix elements not equal to val?
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
virtual Int_t NonZeros() const
Compute the number of elements != 0.0.
virtual Element E2Norm() const
Square of the Euclidian norm, SUM{ m(i,j)^2 }.
Int_t GetRowUpb() const
Definition: TMatrixTBase.h:123
virtual TMatrixTBase< Element > & Apply(const TElementActionT< Element > &action)
Apply action to each matrix element.
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:122
virtual Element Min() const
return minimum matrix element value
virtual Element Max() const
return maximum vector element value
Int_t GetColLwb() const
Definition: TMatrixTBase.h:125
virtual void GetMatrix2Array(Element *data, Option_t *option="") const
Copy matrix data to array .
Int_t GetNoElements() const
Definition: TMatrixTBase.h:128
Bool_t operator<(Element val) const
Are all matrix elements < val?
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.
static Element & NaNValue()
Int_t GetColUpb() const
Definition: TMatrixTBase.h:126
static void DoubleLexSort(Int_t n, Int_t *first, Int_t *second, Element *data)
default kTRUE, when Use array kFALSE
Bool_t operator>(Element val) const
Are all matrix elements > val?
virtual TMatrixTBase< Element > & Sqrt()
Take square root of all elements.
Bool_t IsValid() const
Definition: TMatrixTBase.h:147
Bool_t operator==(Element val) const
Are all matrix elements equal to val?
Int_t GetNcols() const
Definition: TMatrixTBase.h:127
virtual TMatrixTBase< Element > & Abs()
Take an absolute value of a matrix, i.e. apply Abs() to each element.
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.
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.
Bool_t operator<=(Element val) const
Are all matrix elements <= val?
virtual Bool_t IsSymmetric() const
Check whether matrix is symmetric.
virtual TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *option="")
Copy array data to matrix .
Bool_t operator>=(Element val) const
Are all matrix elements >= val?
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 > & NormByDiag(const TVectorT< Element > &v, Option_t *option="D")
option:
Basic string class.
Definition: TString.h:131
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1138
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
TVectorT.
Definition: TVectorT.h:27
double beta(double x, double y)
Calculates the beta function.
const Int_t n
Definition: legend1.C:16
static constexpr double second
static constexpr double m2
Long64_t LocMin(Long64_t n, const T *a)
Return index of array with the minimum element.
Definition: TMath.h:960
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
Long64_t LocMax(Long64_t n, const T *a)
Return index of array with the maximum element.
Definition: TMath.h:988
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
Double_t Log10(Double_t x)
Definition: TMath.h:752
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
Definition: first.py:1
auto * m
Definition: textangle.C:8
static long int sum(long int i)
Definition: Factory.cxx:2258
REAL epsilon
Definition: triangle.c:617