Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
15TMatrixTBase
16
17Template of base class in the linear algebra package.
18
19See the \ref Matrix page for the documentation of the linear algebra package
20
21Matrix properties are stored here, however the data storage is part
22of the derived classes
23*/
24
25#include "TMatrixTBase.h"
26#include "TVectorT.h"
27#include "TBuffer.h"
28#include "TROOT.h"
29#include "TMath.h"
30#include "snprintf.h"
31#include <climits>
32
34
35
36
37////////////////////////////////////////////////////////////////////////////////
38/// Lexical sort on array data using indices first and second
39
40template<class Element>
42{
43 const int incs[] = {1,5,19,41,109,209,505,929,2161,3905,8929,16001,INT_MAX};
44
45 Int_t kinc = 0;
46 while (incs[kinc] <= n/2)
47 kinc++;
48 kinc -= 1;
49
50 // incs[kinc] is the greatest value in the sequence that is also <= n/2.
51 // If n == {0,1}, kinc == -1 and so no sort will take place.
52
53 for( ; kinc >= 0; kinc--) {
54 const Int_t inc = incs[kinc];
55
56 for (Int_t k = inc; k < n; k++) {
57 const Element tmp = data[k];
58 const Int_t fi = first [k];
59 const Int_t se = second[k];
60 Int_t j;
61 for (j = k; j >= inc; j -= inc) {
62 if ( fi < first[j-inc] || (fi == first[j-inc] && se < second[j-inc]) ) {
63 data [j] = data [j-inc];
64 first [j] = first [j-inc];
65 second[j] = second[j-inc];
66 } else
67 break;
68 }
69 data [j] = tmp;
70 first [j] = fi;
71 second[j] = se;
72 }
73 }
74}
75
76////////////////////////////////////////////////////////////////////////////////
77/// Lexical sort on array data using indices first and second
78
79template<class Element>
82{
83 const int incs[] = {1,5,19,41,109,209,505,929,2161,3905,8929,16001,INT_MAX};
84
85 Int_t kinc = 0;
86 while (incs[kinc] <= n/2)
87 kinc++;
88 kinc -= 1;
89
90 // incs[kinc] is the greatest value in the sequence that is also less
91 // than n/2.
92
93 for( ; kinc >= 0; kinc--) {
94 const Int_t inc = incs[kinc];
95
96 if ( !swapFirst && !swapSecond ) {
97 for (Int_t k = inc; k < n; k++) {
98 // loop over all subarrays defined by the current increment
99 const Int_t ktemp = index[k];
100 const Int_t fi = first [ktemp];
101 const Int_t se = second[ktemp];
102 // Insert element k into the sorted subarray
104 for (j = k; j >= inc; j -= inc) {
105 // Loop over the elements in the current subarray
106 if (fi < first[index[j-inc]] || (fi == first[index[j-inc]] && se < second[index[j-inc]])) {
107 // Swap elements j and j - inc, implicitly use the fact
108 // that ktemp hold element j to avoid having to assign to
109 // element j-inc
110 index[j] = index[j-inc];
111 } else {
112 // There are no more elements in this sorted subarray which
113 // are less than element j
114 break;
115 }
116 } // End loop over the elements in the current subarray
117 // Move index[j] out of temporary storage
118 index[j] = ktemp;
119 // The element has been inserted into the subarray.
120 } // End loop over all subarrays defined by the current increment
121 } else if ( swapSecond && !swapFirst ) {
122 for (Int_t k = inc; k < n; k++) {
123 const Int_t ktemp = index[k];
124 const Int_t fi = first [ktemp];
125 const Int_t se = second[k];
126 Int_t j;
127 for (j = k; j >= inc; j -= inc) {
128 if (fi < first[index[j-inc]] || (fi == first[index[j-inc]] && se < second[j-inc])) {
129 index [j] = index[j-inc];
130 second[j] = second[j-inc];
131 } else {
132 break;
133 }
134 }
135 index[j] = ktemp;
136 second[j] = se;
137 }
138 } else if (swapFirst && !swapSecond) {
139 for (Int_t k = inc; k < n; k++ ) {
140 const Int_t ktemp = index[k];
141 const Int_t fi = first[k];
142 const Int_t se = second[ktemp];
143 Int_t j;
144 for (j = k; j >= inc; j -= inc) {
145 if ( fi < first[j-inc] || (fi == first[j-inc] && se < second[ index[j-inc]])) {
146 index[j] = index[j-inc];
147 first[j] = first[j-inc];
148 } else {
149 break;
150 }
151 }
152 index[j] = ktemp;
153 first[j] = fi;
155 } else { // Swap both
156 for (Int_t k = inc; k < n; k++ ) {
157 const Int_t ktemp = index[k];
158 const Int_t fi = first [k];
159 const Int_t se = second[k];
160 Int_t j;
161 for (j = k; j >= inc; j -= inc) {
162 if ( fi < first[j-inc] || (fi == first[j-inc] && se < second[j-inc])) {
163 index [j] = index [j-inc];
164 first [j] = first [j-inc];
165 second[j] = second[j-inc];
166 } else {
167 break;
170 index[j] = ktemp;
171 first[j] = fi;
172 second[j] = se;
176}
177
178////////////////////////////////////////////////////////////////////////////////
179/// Copy array data to matrix . It is assumed that array is of size >= fNelems
180/// (=)))) fNrows*fNcols
181/// option indicates how the data is stored in the array:
182/// option =
183/// - 'F' : column major (Fortran) m[i][j] = array[i+j*fNrows]
184/// - else : row major (C) m[i][j] = array[i*fNcols+j] (default)
185
186template<class Element>
188{
189 R__ASSERT(IsValid());
192 opt.ToUpper();
194 Element *elem = GetMatrixArray();
195 if (opt.Contains("F")) {
196 for (Int_t irow = 0; irow < fNrows; irow++) {
197 const Int_t off1 = irow*fNcols;
198 Int_t off2 = 0;
199 for (Int_t icol = 0; icol < fNcols; icol++) {
200 elem[off1+icol] = data[off2+irow];
201 off2 += fNrows;
203 }
205 else
206 memcpy(elem,data,fNelems*sizeof(Element));
207
208 return *this;
209}
210
211////////////////////////////////////////////////////////////////////////////////
212/// Check whether matrix is symmetric
213
214template<class Element>
216{
217 R__ASSERT(IsValid());
218
219 if ((fNrows != fNcols) || (fRowLwb != fColLwb))
220 return kFALSE;
221
222 const Element * const elem = GetMatrixArray();
223 for (Int_t irow = 0; irow < fNrows; irow++) {
224 const Int_t rowOff = irow*fNcols;
225 Int_t colOff = 0;
226 for (Int_t icol = 0; icol < irow; icol++) {
227 if (elem[rowOff+icol] != elem[colOff+irow])
228 return kFALSE;
229 colOff += fNrows;
230 }
231 }
232 return kTRUE;
233}
234
235////////////////////////////////////////////////////////////////////////////////
236/// Copy matrix data to array . It is assumed that array is of size >= fNelems
237/// (=)))) fNrows*fNcols
238/// option indicates how the data is stored in the array:
239/// option =
240/// - 'F' : column major (Fortran) array[i+j*fNrows] = m[i][j]
241/// - else : row major (C) array[i*fNcols+j] = m[i][j] (default)
242
243template<class Element>
245{
246 R__ASSERT(IsValid());
247
248 TString opt = option;
249 opt.ToUpper();
250
251 const Element * const elem = GetMatrixArray();
252 if (opt.Contains("F")) {
253 for (Int_t irow = 0; irow < fNrows; irow++) {
254 const Int_t off1 = irow*fNcols;
255 Int_t off2 = 0;
256 for (Int_t icol = 0; icol < fNcols; icol++) {
257 data[off2+irow] = elem[off1+icol];
258 off2 += fNrows;
259 }
260 }
261 }
262 else
263 memcpy(data,elem,fNelems*sizeof(Element));
264}
265
266////////////////////////////////////////////////////////////////////////////////
267/// Copy n elements from array v to row rown starting at column coln
268
269template<class Element>
271{
272 const Int_t arown = rown-fRowLwb;
273 const Int_t acoln = coln-fColLwb;
274 const Int_t nr = (n > 0) ? n : fNcols;
275
276 if (gMatrixCheck) {
277 if (arown >= fNrows || arown < 0) {
278 Error("InsertRow","row %d out of matrix range",rown);
279 return *this;
280 }
281
282 if (acoln >= fNcols || acoln < 0) {
283 Error("InsertRow","column %d out of matrix range",coln);
284 return *this;
285 }
286
287 if (acoln+nr > fNcols || nr < 0) {
288 Error("InsertRow","row length %d out of range",nr);
289 return *this;
290 }
291 }
292
293 const Int_t off = arown*fNcols+acoln;
294 Element * const elem = GetMatrixArray()+off;
295 memcpy(elem,v,nr*sizeof(Element));
296
297 return *this;
298}
299
300////////////////////////////////////////////////////////////////////////////////
301/// Store in array v, n matrix elements of row rown starting at column coln
302
303template<class Element>
305{
306 const Int_t arown = rown-fRowLwb;
307 const Int_t acoln = coln-fColLwb;
308 const Int_t nr = (n > 0) ? n : fNcols;
309
310 if (gMatrixCheck) {
311 if (arown >= fNrows || arown < 0) {
312 Error("ExtractRow","row %d out of matrix range",rown);
313 return;
314 }
315
316 if (acoln >= fNcols || acoln < 0) {
317 Error("ExtractRow","column %d out of matrix range",coln);
318 return;
319 }
320
321 if (acoln+n >= fNcols || nr < 0) {
322 Error("ExtractRow","row length %d out of range",nr);
323 return;
324 }
325 }
326
327 const Int_t off = arown*fNcols+acoln;
328 const Element * const elem = GetMatrixArray()+off;
329 memcpy(v,elem,nr*sizeof(Element));
330}
331
332////////////////////////////////////////////////////////////////////////////////
333/// Shift the row index by adding row_shift and the column index by adding
334/// col_shift, respectively. So [rowLwb..rowUpb][colLwb..colUpb] becomes
335/// [rowLwb+row_shift..rowUpb+row_shift][colLwb+col_shift..colUpb+col_shift]
336
337template<class Element>
339{
340 fRowLwb += row_shift;
341 fColLwb += col_shift;
342
343 return *this;
344}
345
346////////////////////////////////////////////////////////////////////////////////
347/// Set matrix elements to zero
348
349template<class Element>
351{
352 R__ASSERT(IsValid());
353 memset(this->GetMatrixArray(),0,fNelems*sizeof(Element));
354
355 return *this;
356}
357
358////////////////////////////////////////////////////////////////////////////////
359/// Take an absolute value of a matrix, i.e. apply Abs() to each element.
360
361template<class Element>
363{
364 R__ASSERT(IsValid());
365
366 Element *ep = this->GetMatrixArray();
367 const Element * const fp = ep+fNelems;
368 while (ep < fp) {
369 *ep = TMath::Abs(*ep);
370 ep++;
371 }
372
373 return *this;
374}
375
376////////////////////////////////////////////////////////////////////////////////
377/// Square each element of the matrix.
378
379template<class Element>
381{
382 R__ASSERT(IsValid());
383
384 Element *ep = this->GetMatrixArray();
385 const Element * const fp = ep+fNelems;
386 while (ep < fp) {
387 *ep = (*ep) * (*ep);
388 ep++;
389 }
390
391 return *this;
392}
393
394////////////////////////////////////////////////////////////////////////////////
395/// Take square root of all elements.
396
397template<class Element>
399{
400 R__ASSERT(IsValid());
401
402 Element *ep = this->GetMatrixArray();
403 const Element * const fp = ep+fNelems;
404 while (ep < fp) {
405 *ep = TMath::Sqrt(*ep);
406 ep++;
407 }
408
409 return *this;
410}
411
412////////////////////////////////////////////////////////////////////////////////
413/// Make a unit matrix (matrix need not be a square one).
414
415template<class Element>
417{
418 R__ASSERT(IsValid());
419
420 Element *ep = this->GetMatrixArray();
421 memset(ep,0,fNelems*sizeof(Element));
422 for (Int_t i = fRowLwb; i <= fRowLwb+fNrows-1; i++)
423 for (Int_t j = fColLwb; j <= fColLwb+fNcols-1; j++)
424 *ep++ = (i==j ? 1.0 : 0.0);
425
426 return *this;
427}
428
429////////////////////////////////////////////////////////////////////////////////
430/// option:
431/// - "D" : b(i,j) = a(i,j)/sqrt(abs*(v(i)*v(j))) (default)
432/// - else : b(i,j) = a(i,j)*sqrt(abs*(v(i)*v(j))) (default)
433
434template<class Element>
436{
437 R__ASSERT(IsValid());
438 R__ASSERT(v.IsValid());
439
440 if (gMatrixCheck) {
441 const Int_t nMax = TMath::Max(fNrows,fNcols);
442 if (v.GetNoElements() < nMax) {
443 Error("NormByDiag","vector shorter than matrix diagonal");
444 return *this;
445 }
446 }
447
448 TString opt(option);
449 opt.ToUpper();
450 const Int_t divide = (opt.Contains("D")) ? 1 : 0;
451
452 const Element *pV = v.GetMatrixArray();
453 Element *mp = this->GetMatrixArray();
454
455 if (divide) {
456 for (Int_t irow = 0; irow < fNrows; irow++) {
457 if (pV[irow] != 0.0) {
458 for (Int_t icol = 0; icol < fNcols; icol++) {
459 if (pV[icol] != 0.0) {
460 const Element val = TMath::Sqrt(TMath::Abs(pV[irow]*pV[icol]));
461 *mp++ /= val;
462 } else {
463 Error("NormbyDiag","vector element %d is zero",icol);
464 mp++;
465 }
466 }
467 } else {
468 Error("NormbyDiag","vector element %d is zero",irow);
469 mp += fNcols;
470 }
471 }
472 } else {
473 for (Int_t irow = 0; irow < fNrows; irow++) {
474 for (Int_t icol = 0; icol < fNcols; icol++) {
475 const Element val = TMath::Sqrt(TMath::Abs(pV[irow]*pV[icol]));
476 *mp++ *= val;
477 }
478 }
479 }
480
481 return *this;
482}
483
484////////////////////////////////////////////////////////////////////////////////
485/// Row matrix norm, MAX{ SUM{ |M(i,j)|, over j}, over i}.
486/// The norm is induced by the infinity vector norm.
487
488template<class Element>
490{
491 R__ASSERT(IsValid());
492
493 const Element * ep = GetMatrixArray();
494 const Element * const fp = ep+fNelems;
495 Element norm = 0;
496
497 // Scan the matrix row-after-row
498 while (ep < fp) {
499 Element sum = 0;
500 // Scan a row to compute the sum
501 for (Int_t j = 0; j < fNcols; j++)
502 sum += TMath::Abs(*ep++);
504 }
505
506 R__ASSERT(ep == fp);
507
508 return norm;
509}
510
511////////////////////////////////////////////////////////////////////////////////
512/// Column matrix norm, MAX{ SUM{ |M(i,j)|, over i}, over j}.
513/// The norm is induced by the 1 vector norm.
514
515template<class Element>
517{
518 R__ASSERT(IsValid());
519
520 const Element * ep = GetMatrixArray();
521 const Element * const fp = ep+fNcols;
522 Element norm = 0;
523
524 // Scan the matrix col-after-col
525 while (ep < fp) {
526 Element sum = 0;
527 // Scan a col to compute the sum
528 for (Int_t i = 0; i < fNrows; i++,ep += fNcols)
529 sum += TMath::Abs(*ep);
530 ep -= fNelems-1; // Point ep to the beginning of the next col
532 }
533
534 R__ASSERT(ep == fp);
535
536 return norm;
537}
538
539////////////////////////////////////////////////////////////////////////////////
540/// Square of the Euclidean norm, SUM{ m(i,j)^2 }.
541
542template<class Element>
544{
545 R__ASSERT(IsValid());
546
547 const Element * ep = GetMatrixArray();
548 const Element * const fp = ep+fNelems;
549 Element sum = 0;
550
551 for ( ; ep < fp; ep++)
552 sum += (*ep) * (*ep);
553
554 return sum;
555}
556
557////////////////////////////////////////////////////////////////////////////////
558/// Compute the number of elements != 0.0
559
560template<class Element>
562{
563 R__ASSERT(IsValid());
564
565 Int_t nr_nonzeros = 0;
566 const Element *ep = this->GetMatrixArray();
567 const Element * const fp = ep+fNelems;
568 while (ep < fp)
569 if (*ep++ != 0.0) nr_nonzeros++;
570
571 return nr_nonzeros;
572}
573
574////////////////////////////////////////////////////////////////////////////////
575/// Compute sum of elements
576
577template<class Element>
579{
580 R__ASSERT(IsValid());
581
582 Element sum = 0.0;
583 const Element *ep = this->GetMatrixArray();
584 const Element * const fp = ep+fNelems;
585 while (ep < fp)
586 sum += *ep++;
587
588 return sum;
589}
590
591////////////////////////////////////////////////////////////////////////////////
592/// return minimum matrix element value
593
594template<class Element>
596{
597 R__ASSERT(IsValid());
598
599 const Element * const ep = this->GetMatrixArray();
600 const Int_t index = TMath::LocMin(fNelems,ep);
601 return ep[index];
602}
603
604////////////////////////////////////////////////////////////////////////////////
605/// return maximum vector element value
606
607template<class Element>
609{
610 R__ASSERT(IsValid());
611
612 const Element * const ep = this->GetMatrixArray();
613 const Int_t index = TMath::LocMax(fNelems,ep);
614 return ep[index];
615}
616
617////////////////////////////////////////////////////////////////////////////////
618/// Draw this matrix
619/// The histogram is named "TMatrixT" by default and no title
620
621template<class Element>
623{
624 gROOT->ProcessLine(Form("THistPainter::PaintSpecialObjects((TObject*)0x%zx,\"%s\");",
625 (size_t)this, option));
626}
627
628////////////////////////////////////////////////////////////////////////////////
629/// Print the matrix as a table of elements.
630/// By default the format "%11.4g" is used to print one element.
631/// One can specify an alternative format with eg
632/// option ="f= %6.2f "
633
634template<class Element>
636{
637 if (!IsValid()) {
638 Error("Print","Matrix is invalid");
639 return;
640 }
641
642 //build format
643 const char *format = "%11.4g ";
644 if (option) {
645 const char *f = strstr(option,"f=");
646 if (f) format = f+2;
647 }
648 char topbar[100];
649 snprintf(topbar,100,format,123.456789);
650 Int_t nch = strlen(topbar)+1;
651 if (nch > 18) nch = 18;
652 char ftopbar[20];
653 for (Int_t i = 0; i < nch; i++) ftopbar[i] = ' ';
654 Int_t nk = 1 + Int_t(TMath::Log10(fNcols));
655 snprintf(ftopbar+nch/2,20-nch/2,"%s%dd","%",nk);
657 for (Int_t i = nch2; i < nch; i++) ftopbar[i] = ' ';
658 ftopbar[nch] = '|';
659 ftopbar[nch+1] = 0;
660
661 printf("\n%dx%d matrix is as follows",fNrows,fNcols);
662
664 if (nch <= 8) cols_per_sheet =10;
665 const Int_t ncols = fNcols;
666 const Int_t nrows = fNrows;
667 const Int_t collwb = fColLwb;
668 const Int_t rowlwb = fRowLwb;
669 nk = 5+nch*TMath::Min(cols_per_sheet,fNcols);
670 for (Int_t i = 0; i < nk; i++) topbar[i] = '-';
671 topbar[nk] = 0;
673 printf("\n\n |");
676 printf("\n%s\n",topbar);
677 if (fNelems <= 0) continue;
678 for (Int_t i = 1; i <= nrows; i++) {
679 printf("%4d |",i+rowlwb-1);
681 printf(format,(*this)(i+rowlwb-1,j+collwb-1));
682 printf("\n");
683 }
684 }
685 printf("\n");
686}
687
688////////////////////////////////////////////////////////////////////////////////
689/// Are all matrix elements equal to val?
690
691template<class Element>
693{
694 R__ASSERT(IsValid());
695
696 if (val == 0. && fNelems == 0)
697 return kTRUE;
698
699 const Element * ep = GetMatrixArray();
700 const Element * const fp = ep+fNelems;
701 for (; ep < fp; ep++)
702 if (!(*ep == val))
703 return kFALSE;
704
705 return kTRUE;
706}
707
708////////////////////////////////////////////////////////////////////////////////
709/// Are all matrix elements not equal to val?
710
711template<class Element>
713{
714 R__ASSERT(IsValid());
715
716 if (val == 0. && fNelems == 0)
717 return kFALSE;
718
719 const Element * ep = GetMatrixArray();
720 const Element * const fp = ep+fNelems;
721 for (; ep < fp; ep++)
722 if (!(*ep != val))
723 return kFALSE;
724
725 return kTRUE;
726}
727
728////////////////////////////////////////////////////////////////////////////////
729/// Are all matrix elements < val?
730
731template<class Element>
733{
734 R__ASSERT(IsValid());
735
736 const Element * ep = GetMatrixArray();
737 const Element * const fp = ep+fNelems;
738 for (; ep < fp; ep++)
739 if (!(*ep < val))
740 return kFALSE;
741
742 return kTRUE;
743}
744
745////////////////////////////////////////////////////////////////////////////////
746/// Are all matrix elements <= val?
747
748template<class Element>
750{
751 R__ASSERT(IsValid());
752
753 const Element * ep = GetMatrixArray();
754 const Element * const fp = ep+fNelems;
755 for (; ep < fp; ep++)
756 if (!(*ep <= val))
757 return kFALSE;
758
759 return kTRUE;
760}
761
762////////////////////////////////////////////////////////////////////////////////
763/// Are all matrix elements > val?
764
765template<class Element>
767{
768 R__ASSERT(IsValid());
769
770 const Element * ep = GetMatrixArray();
771 const Element * const fp = ep+fNelems;
772 for (; ep < fp; ep++)
773 if (!(*ep > val))
774 return kFALSE;
775
776 return kTRUE;
777}
778
779////////////////////////////////////////////////////////////////////////////////
780/// Are all matrix elements >= val?
781
782template<class Element>
784{
785 R__ASSERT(IsValid());
786
787 const Element * ep = GetMatrixArray();
788 const Element * const fp = ep+fNelems;
789 for (; ep < fp; ep++)
790 if (!(*ep >= val))
791 return kFALSE;
792
793 return kTRUE;
794}
795
796////////////////////////////////////////////////////////////////////////////////
797/// Apply action to each matrix element
798
799template<class Element>
801{
802 R__ASSERT(IsValid());
803
804 Element *ep = this->GetMatrixArray();
805 const Element * const ep_last = ep+fNelems;
806 while (ep < ep_last)
807 action.Operation(*ep++);
808
809 return *this;
810}
811
812////////////////////////////////////////////////////////////////////////////////
813/// Apply action to each element of the matrix. To action the location
814/// of the current element is passed.
815
816template<class Element>
818{
819 R__ASSERT(IsValid());
820
821 Element *ep = this->GetMatrixArray();
822 for (action.fI = fRowLwb; action.fI < fRowLwb+fNrows; action.fI++)
823 for (action.fJ = fColLwb; action.fJ < fColLwb+fNcols; action.fJ++)
824 action.Operation(*ep++);
825
826 R__ASSERT(ep == this->GetMatrixArray()+fNelems);
827
828 return *this;
829}
830
831////////////////////////////////////////////////////////////////////////////////
832/// Randomize matrix element values
833
834template<class Element>
836{
837 R__ASSERT(IsValid());
838
839 const Element scale = beta-alpha;
840 const Element shift = alpha/scale;
841
842 Element * ep = GetMatrixArray();
843 const Element * const fp = ep+fNelems;
844 while (ep < fp)
845 *ep++ = scale*(Drand(seed)+shift);
846
847 return *this;
848}
849
850////////////////////////////////////////////////////////////////////////////////
851/// Check to see if two matrices are identical.
852
853template<class Element>
855{
856 if (!AreCompatible(m1,m2)) return kFALSE;
857 return (memcmp(m1.GetMatrixArray(),m2.GetMatrixArray(),
858 m1.GetNoElements()*sizeof(Element)) == 0);
859}
860
861////////////////////////////////////////////////////////////////////////////////
862/// Square of the Euclidean norm of the difference between two matrices.
863
864template<class Element>
866{
867 if (gMatrixCheck && !AreCompatible(m1,m2)) {
868 ::Error("E2Norm","matrices not compatible");
869 return -1.0;
870 }
871
872 const Element * mp1 = m1.GetMatrixArray();
873 const Element * mp2 = m2.GetMatrixArray();
874 const Element * const fmp1 = mp1+m1.GetNoElements();
875
876 Element sum = 0.0;
877 for (; mp1 < fmp1; mp1++, mp2++)
878 sum += (*mp1 - *mp2)*(*mp1 - *mp2);
879
880 return sum;
881}
882
883////////////////////////////////////////////////////////////////////////////////
884/// Check that matrice sm1 and m2 areboth valid and have identical shapes .
885
886template<class Element1,class Element2>
888{
889 if (!m1.IsValid()) {
890 if (verbose)
891 ::Error("AreCompatible", "matrix 1 not valid");
892 return kFALSE;
893 }
894 if (!m2.IsValid()) {
895 if (verbose)
896 ::Error("AreCompatible", "matrix 2 not valid");
897 return kFALSE;
898 }
899
900 if (m1.GetNrows() != m2.GetNrows() || m1.GetNcols() != m2.GetNcols() ||
901 m1.GetRowLwb() != m2.GetRowLwb() || m1.GetColLwb() != m2.GetColLwb()) {
902 if (verbose)
903 ::Error("AreCompatible", "matrices 1 and 2 not compatible");
904 return kFALSE;
905 }
906
907 return kTRUE;
908}
909
910////////////////////////////////////////////////////////////////////////////////
911/// Compare two matrices and print out the result of the comparison.
912
913template<class Element>
915{
916 if (!AreCompatible(m1,m2)) {
917 Error("Compare(const TMatrixTBase<Element> &,const TMatrixTBase<Element> &)","matrices are incompatible");
918 return;
919 }
920
921 printf("\n\nComparison of two TMatrices:\n");
922
923 Element norm1 = 0; // Norm of the Matrices
924 Element norm2 = 0;
925 Element ndiff = 0; // Norm of the difference
926 Int_t imax = 0; // For the elements that differ most
927 Int_t jmax = 0;
928 Element difmax = -1;
929
930 for (Int_t i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++) {
931 for (Int_t j = m1.GetColLwb(); j < m1.GetColUpb(); j++) {
932 const Element mv1 = m1(i,j);
933 const Element mv2 = m2(i,j);
934 const Element diff = TMath::Abs(mv1-mv2);
935
936 if (diff > difmax) {
937 difmax = diff;
938 imax = i;
939 jmax = j;
940 }
941 norm1 += TMath::Abs(mv1);
942 norm2 += TMath::Abs(mv2);
944 }
945 }
946
947 printf("\nMaximal discrepancy \t\t%g", difmax);
948 printf("\n occurred at the point\t\t(%d,%d)",imax,jmax);
949 const Element mv1 = m1(imax,jmax);
950 const Element mv2 = m2(imax,jmax);
951 printf("\n Matrix 1 element is \t\t%g", mv1);
952 printf("\n Matrix 2 element is \t\t%g", mv2);
953 printf("\n Absolute error v2[i]-v1[i]\t\t%g", mv2-mv1);
954 printf("\n Relative error\t\t\t\t%g\n",
955 (mv2-mv1)/TMath::Max(TMath::Abs(mv2+mv1)/2,(Element)1e-7));
956
957 printf("\n||Matrix 1|| \t\t\t%g", norm1);
958 printf("\n||Matrix 2|| \t\t\t%g", norm2);
959 printf("\n||Matrix1-Matrix2||\t\t\t\t%g", ndiff);
960 printf("\n||Matrix1-Matrix2||/sqrt(||Matrix1|| ||Matrix2||)\t%g\n\n",
962}
963
964////////////////////////////////////////////////////////////////////////////////
965/// Validate that all elements of matrix have value val within maxDevAllow.
966
967template<class Element>
969{
970 R__ASSERT(m.IsValid());
971
972 if (m == 0)
973 return kTRUE;
974
975 Int_t imax = 0;
976 Int_t jmax = 0;
977 Element maxDevObs = 0;
978
979 if (TMath::Abs(maxDevAllow) <= 0.0)
980 maxDevAllow = std::numeric_limits<Element>::epsilon();
981
982 for (Int_t i = m.GetRowLwb(); i <= m.GetRowUpb(); i++) {
983 for (Int_t j = m.GetColLwb(); j <= m.GetColUpb(); j++) {
984 const Element dev = TMath::Abs(m(i,j)-val);
985 if (dev > maxDevObs) {
986 imax = i;
987 jmax = j;
988 maxDevObs = dev;
989 }
990 }
991 }
992
993 if (maxDevObs == 0)
994 return kTRUE;
995
996 if (verbose) {
997 printf("Largest dev for (%d,%d); dev = |%g - %g| = %g\n",imax,jmax,m(imax,jmax),val,maxDevObs);
999 Error("VerifyElementValue","Deviation > %g\n",maxDevAllow);
1000 }
1001
1003 return kFALSE;
1004 return kTRUE;
1005}
1006
1007////////////////////////////////////////////////////////////////////////////////
1008/// Verify that elements of the two matrices are equal within MaxDevAllow .
1009
1010template<class Element>
1012 Element maxDevAllow)
1013{
1014 if (!AreCompatible(m1,m2,verbose))
1015 return kFALSE;
1016
1017 if (m1 == 0 && m2 == 0)
1018 return kTRUE;
1019
1020 Int_t imax = 0;
1021 Int_t jmax = 0;
1022 Element maxDevObs = 0;
1023
1024 if (TMath::Abs(maxDevAllow) <= 0.0)
1025 maxDevAllow = std::numeric_limits<Element>::epsilon();
1026
1027 for (Int_t i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++) {
1028 for (Int_t j = m1.GetColLwb(); j <= m1.GetColUpb(); j++) {
1029 const Element dev = TMath::Abs(m1(i,j)-m2(i,j));
1030 if (dev > maxDevObs) {
1031 imax = i;
1032 jmax = j;
1033 maxDevObs = dev;
1034 }
1035 }
1036 }
1037
1038 if (maxDevObs == 0)
1039 return kTRUE;
1040
1041 if (verbose) {
1042 printf("Largest dev for (%d,%d); dev = |%g - %g| = %g\n",
1043 imax,jmax,m1(imax,jmax),m2(imax,jmax),maxDevObs);
1044 if (maxDevObs > maxDevAllow)
1045 Error("VerifyMatrixValue","Deviation > %g\n",maxDevAllow);
1046 }
1047
1048 if (maxDevObs > maxDevAllow)
1049 return kFALSE;
1050 return kTRUE;
1051}
1052
1053////////////////////////////////////////////////////////////////////////////////
1054/// Stream an object of class TMatrixTBase<Element>.
1055
1056template<class Element>
1058{
1059 if (R__b.IsReading()) {
1060 UInt_t R__s, R__c;
1061 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
1062 if (R__v > 1) {
1063 R__b.ReadClassBuffer(TMatrixTBase<Element>::Class(),this,R__v,R__s,R__c);
1064 } else {
1065 Error("TMatrixTBase<Element>::Streamer","Unknown version number: %d",R__v);
1066 R__ASSERT(0);
1067 }
1068 if (R__v < 4) MakeValid();
1069 } else {
1070 R__b.WriteClassBuffer(TMatrixTBase<Element>::Class(),this);
1071 }
1072}
1073
1074// trick to return a reference to nan in operator(i,j_ when i,j are outside of range
1075template<class Element>
1077 static Element gNanValue;
1078};
1079template<>
1080Double_t nan_value_t<Double_t>::gNanValue = std::numeric_limits<Double_t>::quiet_NaN();
1081template<>
1082Float_t nan_value_t<Float_t>::gNanValue = std::numeric_limits<Float_t>::quiet_NaN();
1083
1084template<class Element>
1089
1090
1091template class TMatrixTBase<Float_t>;
1092
1093template Bool_t TMatrixTAutoloadOps::operator== <Float_t>(const TMatrixFBase &m1,const TMatrixFBase &m2);
1094template Float_t TMatrixTAutoloadOps::E2Norm <Float_t>(const TMatrixFBase &m1,const TMatrixFBase &m2);
1095template Bool_t TMatrixTAutoloadOps::AreCompatible<Float_t,Float_t>
1096 (const TMatrixFBase &m1,const TMatrixFBase &m2,Int_t verbose);
1097template Bool_t TMatrixTAutoloadOps::AreCompatible<Float_t,Double_t>
1098 (const TMatrixFBase &m1,const TMatrixDBase &m2,Int_t verbose);
1099template void TMatrixTAutoloadOps::Compare <Float_t>(const TMatrixFBase &m1,const TMatrixFBase &m2);
1100template Bool_t TMatrixTAutoloadOps::VerifyMatrixValue <Float_t>(const TMatrixFBase &m,Float_t val,Int_t verbose,Float_t maxDevAllow);
1101template Bool_t TMatrixTAutoloadOps::VerifyMatrixValue <Float_t>(const TMatrixFBase &m,Float_t val);
1102template Bool_t TMatrixTAutoloadOps::VerifyMatrixIdentity<Float_t>(const TMatrixFBase &m1,const TMatrixFBase &m2,
1103 Int_t verbose,Float_t maxDevAllowN);
1104template Bool_t TMatrixTAutoloadOps::VerifyMatrixIdentity<Float_t>(const TMatrixFBase &m1,const TMatrixFBase &m2);
1105
1106template class TMatrixTBase<Double_t>;
1107
1108template Bool_t TMatrixTAutoloadOps::operator== <Double_t>(const TMatrixDBase &m1,const TMatrixDBase &m2);
1109template Double_t TMatrixTAutoloadOps::E2Norm <Double_t>(const TMatrixDBase &m1,const TMatrixDBase &m2);
1110template Bool_t TMatrixTAutoloadOps::AreCompatible<Double_t,Double_t>
1111 (const TMatrixDBase &m1,const TMatrixDBase &m2,Int_t verbose);
1112template Bool_t TMatrixTAutoloadOps::AreCompatible<Double_t,Float_t>
1113 (const TMatrixDBase &m1,const TMatrixFBase &m2,Int_t verbose);
1114template void TMatrixTAutoloadOps::Compare <Double_t>(const TMatrixDBase &m1,const TMatrixDBase &m2);
1115template Bool_t TMatrixTAutoloadOps::VerifyMatrixValue <Double_t>(const TMatrixDBase &m,Double_t val,Int_t verbose,Double_t maxDevAllow);
1116template Bool_t TMatrixTAutoloadOps::VerifyMatrixValue <Double_t>(const TMatrixDBase &m,Double_t val);
1117template Bool_t TMatrixTAutoloadOps::VerifyMatrixIdentity<Double_t>(const TMatrixDBase &m1,const TMatrixDBase &m2,
1118 Int_t verbose,Double_t maxDevAllow);
1119template Bool_t TMatrixTAutoloadOps::VerifyMatrixIdentity<Double_t>(const TMatrixDBase &m1,const TMatrixDBase &m2);
#define f(i)
Definition RSha256.hxx:104
#define e(i)
Definition RSha256.hxx:103
bool Bool_t
Boolean (0=false, 1=true) (bool)
Definition RtypesCore.h:77
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:59
short Version_t
Class version identifier (short)
Definition RtypesCore.h:79
float Float_t
Float 4 bytes (float)
Definition RtypesCore.h:71
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
const char Option_t
Option string (const char)
Definition RtypesCore.h:80
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
#define R__ASSERT(e)
Checks condition e and reports a fatal error if it's false.
Definition TError.h:125
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:208
Option_t Option_t option
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t format
Int_t gMatrixCheck
R__EXTERN Int_t gMatrixCheck
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
#define gROOT
Definition TROOT.h:411
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2495
#define snprintf
Definition civetweb.c:1579
Buffer base class used for serializing objects.
Definition TBuffer.h:43
TMatrixTBase.
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}.
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 TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
Bool_t operator!=(Element val) const
Are all matrix elements not equal to val?
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
void Print(Option_t *name="") const override
Print the matrix as a table of elements.
virtual Int_t NonZeros() const
Compute the number of elements != 0.0.
virtual Element E2Norm() const
Square of the Euclidean norm, SUM{ m(i,j)^2 }.
virtual TMatrixTBase< Element > & Apply(const TElementActionT< Element > &action)
Apply action to each matrix element.
virtual Element Min() const
return minimum matrix element value
virtual Element Max() const
return maximum vector element value
virtual void GetMatrix2Array(Element *data, Option_t *option="") const
Copy matrix data to array .
Bool_t operator<(Element val) const
Are all matrix elements < val?
void Draw(Option_t *option="") override
Draw this matrix The histogram is named "TMatrixT" by default and no title.
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()
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 operator==(Element val) const
Are all matrix elements equal to val?
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?
void Streamer(TBuffer &) override
Stream an object of class TMatrixTBase<Element>.
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:138
void ToUpper()
Change string to upper case.
Definition TString.cxx:1202
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:640
const Int_t n
Definition legend1.C:16
Long64_t LocMin(Long64_t n, const T *a)
Returns index of array with the minimum element.
Definition TMath.h:993
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:251
Long64_t LocMax(Long64_t n, const T *a)
Returns index of array with the maximum element.
Definition TMath.h:1099
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:199
Double_t Log10(Double_t x)
Returns the common (base-10) logarithm of x.
Definition TMath.h:773
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124
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 .
void Compare(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Compare two matrices and print out the result of the comparison.
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.
Bool_t operator==(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Check to see if two matrices are identical.
Element E2Norm(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Square of the Euclidean norm of the difference between two matrices.
Bool_t AreCompatible(const TMatrixTBase< Element1 > &m1, const TMatrixTBase< Element2 > &m2, Int_t verbose=0)
Check that matrice sm1 and m2 areboth valid and have identical shapes .
static Element gNanValue
TMarker m
Definition textangle.C:8
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2339