Logo ROOT   master
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 
41 template<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 
80 template<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 
187 template<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 
215 template<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 
244 template<class Element>
245 void 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 
270 template<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 
304 template<class Element>
305 void 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 
338 template<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 
350 template<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 
362 template<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 
380 template<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 
398 template<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 
416 template<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 
435 template<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 
489 template<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 
516 template<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 
543 template<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 
561 template<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 
578 template<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 
595 template<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 
608 template<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 
622 template<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 
635 template<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 
692 template<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 
712 template<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 
732 template<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 
749 template<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 
766 template<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 
783 template<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 
800 template<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 
817 template<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 
835 template<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 
854 template<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 
865 template<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 
887 template<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 
914 template<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 
968 template<class Element>
969 Bool_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 
1011 template<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)
1026  maxDevAllow = std::numeric_limits<Element>::epsilon();
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 
1057 template<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
1076 template<class Element>
1077 struct nan_value_t {
1078  static Element gNanValue;
1079 };
1080 template<>
1081 Double_t nan_value_t<Double_t>::gNanValue = std::numeric_limits<Double_t>::quiet_NaN();
1082 template<>
1083 Float_t nan_value_t<Float_t>::gNanValue = std::numeric_limits<Float_t>::quiet_NaN();
1084 
1085 template<class Element>
1087 {
1088  return nan_value_t<Element>::gNanValue;
1089 }
1090 
1091 
1092 template class TMatrixTBase<Float_t>;
1093 
1094 template Bool_t operator== <Float_t>(const TMatrixFBase &m1,const TMatrixFBase &m2);
1095 template Float_t E2Norm <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);
1100 template void Compare <Float_t>(const TMatrixFBase &m1,const TMatrixFBase &m2);
1101 template Bool_t VerifyMatrixValue <Float_t>(const TMatrixFBase &m,Float_t val,Int_t verbose,Float_t maxDevAllow);
1104  Int_t verbose,Float_t maxDevAllowN);
1106 
1107 template class TMatrixTBase<Double_t>;
1108 
1109 template Bool_t operator== <Double_t>(const TMatrixDBase &m1,const TMatrixDBase &m2);
1110 template Double_t E2Norm <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);
1115 template void Compare <Double_t>(const TMatrixDBase &m1,const TMatrixDBase &m2);
1119  Int_t verbose,Double_t maxDevAllow);
Bool_t VerifyMatrixIdentity(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2, Int_t verbose, Element maxDevAllow)
Verify that elements of the two matrices are equal within MaxDevAllow .
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
virtual const Element * GetMatrixArray() const =0
virtual TMatrixTBase< Element > & Sqrt()
Take square root of all elements.
virtual Bool_t IsSymmetric() const
Check whether matrix is symmetric.
Bool_t IsReading() const
Definition: TBuffer.h:85
virtual TMatrixTBase< Element > & Shift(Int_t row_shift, Int_t col_shift)
Shift the row index by adding row_shift and the column index by adding col_shift, respectively...
virtual Int_t NonZeros() const
Compute the number of elements != 0.0.
static long int sum(long int i)
Definition: Factory.cxx:2275
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
template Bool_t AreCompatible< Float_t, Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2, Int_t verbose)
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:122
auto * m
Definition: textangle.C:8
Long64_t LocMax(Long64_t n, const T *a)
Return index of array with the maximum element.
Definition: TMath.h:990
virtual TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *option="")
Copy array data to matrix .
short Version_t
Definition: RtypesCore.h:63
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
Bool_t operator!=(Element val) const
Are all matrix elements not equal to val?
float Float_t
Definition: RtypesCore.h:55
const char Option_t
Definition: RtypesCore.h:64
Int_t GetNcols() const
Definition: TMatrixTBase.h:127
TVectorT.
Definition: TMatrixTBase.h:79
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1138
Buffer base class used for serializing objects.
Definition: TBuffer.h:42
#define R__ASSERT(e)
Definition: TError.h:96
#define gROOT
Definition: TROOT.h:405
Basic string class.
Definition: TString.h:131
#define f(i)
Definition: RSha256.hxx:104
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
virtual TMatrixTBase< Element > & Sqr()
Square each element of the matrix.
int Int_t
Definition: RtypesCore.h:43
bool Bool_t
Definition: RtypesCore.h:61
Bool_t operator>(Element val) const
Are all matrix elements > val?
void Compare(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Compare two matrices and print out the result of the comparison.
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
double beta(double x, double y)
Calculates the beta function.
template Double_t E2Norm< Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2)
template Bool_t VerifyMatrixValue< Float_t >(const TMatrixFBase &m, Float_t val, Int_t verbose, Float_t maxDevAllow)
Int_t GetNoElements() const
Definition: TMatrixTBase.h:128
Bool_t operator<(Element val) const
Are all matrix elements < val?
static void DoubleLexSort(Int_t n, Int_t *first, Int_t *second, Element *data)
default kTRUE, when Use array kFALSE
#define templateClassImp(name)
Definition: Rtypes.h:405
Int_t GetColLwb() const
Definition: TMatrixTBase.h:125
virtual void Operation(Element &element) const =0
Element E2Norm(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Square of the Euclidian norm of the difference between two matrices.
Double_t Log10(Double_t x)
Definition: TMath.h:754
Bool_t operator>=(Element val) const
Are all matrix elements >= val?
template Bool_t AreCompatible< Double_t, Float_t >(const TMatrixDBase &m1, const TMatrixFBase &m2, Int_t verbose)
void Error(const char *location, const char *msgfmt,...)
Bool_t operator==(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Check to see if two matrices are identical.
template Bool_t VerifyMatrixIdentity< Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2, Int_t verbose, Double_t maxDevAllow)
virtual TMatrixTBase< Element > & Apply(const TElementActionT< Element > &action)
Apply action to each matrix element.
virtual TMatrixTBase< Element > & NormByDiag(const TVectorT< Element > &v, Option_t *option="D")
option:
static constexpr double m2
template Bool_t AreCompatible< Double_t, Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2, Int_t verbose)
virtual void ExtractRow(Int_t row, Int_t col, Element *v, Int_t n=-1) const
Store in array v, n matrix elements of row rown starting at column coln.
virtual TMatrixTBase< Element > & InsertRow(Int_t row, Int_t col, const Element *v, Int_t n=-1)
Copy n elements from array v to row rown starting at column coln.
unsigned int UInt_t
Definition: RtypesCore.h:44
char * Form(const char *fmt,...)
virtual void GetMatrix2Array(Element *data, Option_t *option="") const
Copy matrix data to array .
virtual Element Min() const
return minimum matrix element value
virtual TMatrixTBase< Element > & Randomize(Element alpha, Element beta, Double_t &seed)
Randomize matrix element values.
template Bool_t AreCompatible< Float_t, Double_t >(const TMatrixFBase &m1, const TMatrixDBase &m2, Int_t verbose)
Bool_t operator<=(Element val) const
Are all matrix elements <= val?
Int_t GetColUpb() const
Definition: TMatrixTBase.h:126
TMatrixTBase.
REAL epsilon
Definition: triangle.c:617
Int_t GetNrows() const
Definition: TMatrixTBase.h:124
const Bool_t kFALSE
Definition: RtypesCore.h:90
virtual Element RowNorm() const
Row matrix norm, MAX{ SUM{ |M(i,j)|, over j}, over i}.
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
virtual Element E2Norm() const
Square of the Euclidian norm, SUM{ m(i,j)^2 }.
template Bool_t VerifyMatrixIdentity< Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2, Int_t verbose, Float_t maxDevAllowN)
virtual Element Sum() const
Compute sum of elements.
template void Compare< Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2)
double Double_t
Definition: RtypesCore.h:57
template void Compare< Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2)
unsigned long ULong_t
Definition: RtypesCore.h:53
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Bool_t AreCompatible(const TMatrixTBase< Element1 > &m1, const TMatrixTBase< Element2 > &m2, Int_t verbose)
Check that matrice sm1 and m2 areboth valid and have identical shapes .
virtual void Operation(Element &element) const =0
void Print(Option_t *name="") const
Print the matrix as a table of elements.
void Draw(Option_t *option="")
Draw this matrix The histogram is named "TMatrixT" by default and no title.
Int_t gMatrixCheck
static void IndexedLexSort(Int_t n, Int_t *first, Int_t swapFirst, Int_t *second, Int_t swapSecond, Int_t *index)
Lexical sort on array data using indices first and second.
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
Bool_t operator==(Element val) const
Are all matrix elements equal to val?
#define snprintf
Definition: civetweb.c:1540
static Element & NaNValue()
Int_t GetRowUpb() const
Definition: TMatrixTBase.h:123
Definition: first.py:1
virtual Element ColNorm() const
Column matrix norm, MAX{ SUM{ |M(i,j)|, over i}, over j}.
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
Long64_t LocMin(Long64_t n, const T *a)
Return index of array with the minimum element.
Definition: TMath.h:962
const Bool_t kTRUE
Definition: RtypesCore.h:89
template Bool_t VerifyMatrixValue< Double_t >(const TMatrixDBase &m, Double_t val, Int_t verbose, Double_t maxDevAllow)
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
virtual Element Max() const
return maximum vector element value
Bool_t VerifyMatrixValue(const TMatrixTBase< Element > &m, Element val, Int_t verbose, Element maxDevAllow)
Validate that all elements of matrix have value val within maxDevAllow.
template Float_t E2Norm< Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2)
const Int_t n
Definition: legend1.C:16
Bool_t IsValid() const
Definition: TMatrixTBase.h:147
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
virtual TMatrixTBase< Element > & Abs()
Take an absolute value of a matrix, i.e. apply Abs() to each element.
static constexpr double second