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