41template<
class Element>
44 const int incs[] = {1,5,19,41,109,209,505,929,2161,3905,8929,16001,INT_MAX};
47 while (incs[kinc] <=
n/2)
54 for( ; kinc >= 0; kinc--) {
55 const Int_t inc = incs[kinc];
57 for (
Int_t k = inc; k <
n; k++) {
58 const Element tmp = data[k];
60 const Int_t se = second[k];
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];
66 second[j] = second[j-inc];
80template<
class Element>
84 const int incs[] = {1,5,19,41,109,209,505,929,2161,3905,8929,16001,INT_MAX};
87 while (incs[kinc] <=
n/2)
94 for( ; kinc >= 0; kinc--) {
95 const Int_t inc = incs[kinc];
97 if ( !swapFirst && !swapSecond ) {
98 for (
Int_t k = inc; k <
n; k++) {
100 const Int_t ktemp = index[k];
105 for (j = k; j >= inc; j -= inc) {
107 if (fi <
first[index[j-inc]] || (fi ==
first[index[j-inc]] && se < second[index[j-inc]])) {
111 index[j] = index[j-inc];
122 }
else if ( swapSecond && !swapFirst ) {
123 for (
Int_t k = inc; k <
n; k++) {
124 const Int_t ktemp = index[k];
126 const Int_t se = second[k];
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];
139 }
else if (swapFirst && !swapSecond) {
140 for (
Int_t k = inc; k <
n; k++ ) {
141 const Int_t ktemp = index[k];
143 const Int_t se = second[ktemp];
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];
157 for (
Int_t k = inc; k <
n; k++ ) {
160 const Int_t se = second[k];
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];
166 second[j] = second[j-inc];
187template<
class Element>
195 Element *elem = GetMatrixArray();
197 for (
Int_t irow = 0; irow < fNrows; irow++) {
198 const Int_t off1 = irow*fNcols;
200 for (
Int_t icol = 0; icol < fNcols; icol++) {
201 elem[off1+icol] = data[off2+irow];
207 memcpy(elem,data,fNelems*
sizeof(Element));
215template<
class Element>
220 if ((fNrows != fNcols) || (fRowLwb != fColLwb))
223 const Element *
const elem = GetMatrixArray();
224 for (
Int_t irow = 0; irow < fNrows; irow++) {
225 const Int_t rowOff = irow*fNcols;
227 for (
Int_t icol = 0; icol < irow; icol++) {
228 if (elem[rowOff+icol] != elem[colOff+irow])
244template<
class Element>
252 const Element *
const elem = GetMatrixArray();
254 for (
Int_t irow = 0; irow < fNrows; irow++) {
255 const Int_t off1 = irow*fNcols;
257 for (
Int_t icol = 0; icol < fNcols; icol++) {
258 data[off2+irow] = elem[off1+icol];
264 memcpy(data,elem,fNelems*
sizeof(Element));
270template<
class Element>
273 const Int_t arown = rown-fRowLwb;
274 const Int_t acoln = coln-fColLwb;
275 const Int_t nr = (
n > 0) ?
n : fNcols;
278 if (arown >= fNrows || arown < 0) {
279 Error(
"InsertRow",
"row %d out of matrix range",rown);
283 if (acoln >= fNcols || acoln < 0) {
284 Error(
"InsertRow",
"column %d out of matrix range",coln);
288 if (acoln+nr > fNcols || nr < 0) {
289 Error(
"InsertRow",
"row length %d out of range",nr);
294 const Int_t off = arown*fNcols+acoln;
295 Element *
const elem = GetMatrixArray()+off;
296 memcpy(elem,
v,nr*
sizeof(Element));
304template<
class Element>
307 const Int_t arown = rown-fRowLwb;
308 const Int_t acoln = coln-fColLwb;
309 const Int_t nr = (
n > 0) ?
n : fNcols;
312 if (arown >= fNrows || arown < 0) {
313 Error(
"ExtractRow",
"row %d out of matrix range",rown);
317 if (acoln >= fNcols || acoln < 0) {
318 Error(
"ExtractRow",
"column %d out of matrix range",coln);
322 if (acoln+
n >= fNcols || nr < 0) {
323 Error(
"ExtractRow",
"row length %d out of range",nr);
328 const Int_t off = arown*fNcols+acoln;
329 const Element *
const elem = GetMatrixArray()+off;
330 memcpy(
v,elem,nr*
sizeof(Element));
338template<
class Element>
341 fRowLwb += row_shift;
342 fColLwb += col_shift;
350template<
class Element>
354 memset(this->GetMatrixArray(),0,fNelems*
sizeof(Element));
362template<
class Element>
367 Element *ep = this->GetMatrixArray();
368 const Element *
const fp = ep+fNelems;
380template<
class Element>
385 Element *ep = this->GetMatrixArray();
386 const Element *
const fp = ep+fNelems;
398template<
class Element>
403 Element *ep = this->GetMatrixArray();
404 const Element *
const fp = ep+fNelems;
416template<
class Element>
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);
435template<
class Element>
443 if (
v.GetNoElements() < nMax) {
444 Error(
"NormByDiag",
"vector shorter than matrix diagonal");
453 const Element *pV =
v.GetMatrixArray();
454 Element *mp = this->GetMatrixArray();
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) {
464 Error(
"NormbyDiag",
"vector element %d is zero",icol);
469 Error(
"NormbyDiag",
"vector element %d is zero",irow);
474 for (
Int_t irow = 0; irow < fNrows; irow++) {
475 for (
Int_t icol = 0; icol < fNcols; icol++) {
489template<
class Element>
494 const Element * ep = GetMatrixArray();
495 const Element *
const fp = ep+fNelems;
502 for (
Int_t j = 0; j < fNcols; j++)
516template<
class Element>
521 const Element * ep = GetMatrixArray();
522 const Element *
const fp = ep+fNcols;
529 for (
Int_t i = 0; i < fNrows; i++,ep += fNcols)
543template<
class Element>
548 const Element * ep = GetMatrixArray();
549 const Element *
const fp = ep+fNelems;
552 for ( ; ep < fp; ep++)
553 sum += (*ep) * (*ep);
561template<
class Element>
566 Int_t nr_nonzeros = 0;
567 const Element *ep = this->GetMatrixArray();
568 const Element *
const fp = ep+fNelems;
570 if (*ep++ != 0.0) nr_nonzeros++;
578template<
class Element>
584 const Element *ep = this->GetMatrixArray();
585 const Element *
const fp = ep+fNelems;
595template<
class Element>
600 const Element *
const ep = this->GetMatrixArray();
608template<
class Element>
613 const Element *
const ep = this->GetMatrixArray();
622template<
class Element>
625 gROOT->ProcessLine(
Form(
"THistPainter::PaintSpecialObjects((TObject*)0x%zx,\"%s\");",
626 (
size_t)
this, option));
635template<
class Element>
639 Error(
"Print",
"Matrix is invalid");
644 const char *format =
"%11.4g ";
646 const char *
f = strstr(option,
"f=");
650 snprintf(topbar,100,format,123.456789);
651 Int_t nch = strlen(topbar)+1;
652 if (nch > 18) nch = 18;
654 for (
Int_t i = 0; i < nch; i++) ftopbar[i] =
' ';
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] =
' ';
662 printf(
"\n%dx%d matrix is as follows",fNrows,fNcols);
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;
671 for (
Int_t i = 0; i < nk; i++) topbar[i] =
'-';
673 for (
Int_t sheet_counter = 1; sheet_counter <= ncols; sheet_counter += cols_per_sheet) {
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));
692template<
class Element>
697 if (val == 0. && fNelems == 0)
700 const Element * ep = GetMatrixArray();
701 const Element *
const fp = ep+fNelems;
702 for (; ep < fp; ep++)
712template<
class Element>
717 if (val == 0. && fNelems == 0)
720 const Element * ep = GetMatrixArray();
721 const Element *
const fp = ep+fNelems;
722 for (; ep < fp; ep++)
732template<
class Element>
737 const Element * ep = GetMatrixArray();
738 const Element *
const fp = ep+fNelems;
739 for (; ep < fp; ep++)
749template<
class Element>
754 const Element * ep = GetMatrixArray();
755 const Element *
const fp = ep+fNelems;
756 for (; ep < fp; ep++)
766template<
class Element>
771 const Element * ep = GetMatrixArray();
772 const Element *
const fp = ep+fNelems;
773 for (; ep < fp; ep++)
783template<
class Element>
788 const Element * ep = GetMatrixArray();
789 const Element *
const fp = ep+fNelems;
790 for (; ep < fp; ep++)
800template<
class Element>
805 Element *ep = this->GetMatrixArray();
806 const Element *
const ep_last = ep+fNelems;
817template<
class Element>
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++)
827 R__ASSERT(ep == this->GetMatrixArray()+fNelems);
835template<
class Element>
840 const Element scale = beta-alpha;
841 const Element shift = alpha/scale;
843 Element * ep = GetMatrixArray();
844 const Element *
const fp = ep+fNelems;
846 *ep++ = scale*(
Drand(seed)+shift);
854template<
class Element>
865template<
class Element>
869 ::Error(
"E2Norm",
"matrices not compatible");
878 for (; mp1 < fmp1; mp1++, mp2++)
879 sum += (*mp1 - *mp2)*(*mp1 - *mp2);
887template<
class Element1,
class Element2>
892 ::Error(
"AreCompatible",
"matrix 1 not valid");
897 ::Error(
"AreCompatible",
"matrix 2 not valid");
904 ::Error(
"AreCompatible",
"matrices 1 and 2 not compatible");
914template<
class Element>
918 Error(
"Compare(const TMatrixTBase<Element> &,const TMatrixTBase<Element> &)",
"matrices are incompatible");
922 printf(
"\n\nComparison of two TMatrices:\n");
933 const Element mv1 = m1(i,j);
934 const Element mv2 = m2(i,j);
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",
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",
968template<
class Element>
978 Element maxDevObs = 0;
981 maxDevAllow = std::numeric_limits<Element>::epsilon();
983 for (
Int_t i =
m.GetRowLwb(); i <=
m.GetRowUpb(); i++) {
984 for (
Int_t j =
m.GetColLwb(); j <=
m.GetColUpb(); j++) {
986 if (dev > maxDevObs) {
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);
1003 if(maxDevObs > maxDevAllow)
1011template<
class Element>
1013 Element maxDevAllow)
1018 if (m1 == 0 && m2 == 0)
1023 Element maxDevObs = 0;
1026 maxDevAllow = std::numeric_limits<Element>::epsilon();
1030 const Element dev =
TMath::Abs(m1(i,j)-m2(i,j));
1031 if (dev > maxDevObs) {
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);
1049 if (maxDevObs > maxDevAllow)
1057template<
class Element>
1066 Error(
"TMatrixTBase<Element>::Streamer",
"Unknown version number: %d",R__v);
1069 if (R__v < 4) MakeValid();
1076template<
class Element>
1085template<
class Element>
1096template Bool_t TMatrixTAutoloadOps::AreCompatible<Float_t,Float_t>
1098template Bool_t TMatrixTAutoloadOps::AreCompatible<Float_t,Double_t>
1111template Bool_t TMatrixTAutoloadOps::AreCompatible<Double_t,Double_t>
1113template Bool_t TMatrixTAutoloadOps::AreCompatible<Double_t,Float_t>
#define templateClassImp(name)
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
char * Form(const char *fmt,...)
Buffer base class used for serializing objects.
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
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
virtual void Operation(Element &element) const =0
virtual void Operation(Element &element) const =0
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
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 }.
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 .
Int_t GetNoElements() const
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()
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?
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:
void ToUpper()
Change string to upper case.
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Long64_t LocMin(Long64_t n, const T *a)
Return index of array with the minimum element.
Short_t Max(Short_t a, Short_t b)
Long64_t LocMax(Long64_t n, const T *a)
Return index of array with the maximum element.
Double_t Sqrt(Double_t x)
Short_t Min(Short_t a, Short_t b)
Double_t Log10(Double_t x)
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 Euclidian 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 uint64_t sum(uint64_t i)