237template<
class Element>
240 const int incs[] = {1,5,19,41,109,209,505,929,2161,3905,8929,16001,INT_MAX};
243 while (incs[kinc] <=
n/2)
250 for( ; kinc >= 0; kinc--) {
251 const Int_t inc = incs[kinc];
253 for (
Int_t k = inc; k <
n; k++) {
254 const Element tmp =
data[k];
258 for (j = k; j >= inc; j -= inc) {
276template<
class Element>
280 const int incs[] = {1,5,19,41,109,209,505,929,2161,3905,8929,16001,INT_MAX};
283 while (incs[kinc] <=
n/2)
290 for( ; kinc >= 0; kinc--) {
291 const Int_t inc = incs[kinc];
293 if ( !swapFirst && !swapSecond ) {
294 for (
Int_t k = inc; k <
n; k++) {
296 const Int_t ktemp = index[k];
301 for (j = k; j >= inc; j -= inc) {
303 if (fi <
first[index[j-inc]] || (fi ==
first[index[j-inc]] && se <
second[index[j-inc]])) {
307 index[j] = index[j-inc];
318 }
else if ( swapSecond && !swapFirst ) {
319 for (
Int_t k = inc; k <
n; k++) {
320 const Int_t ktemp = index[k];
324 for (j = k; j >= inc; j -= inc) {
325 if (fi <
first[index[j-inc]] || (fi ==
first[index[j-inc]] && se <
second[j-inc])) {
326 index [j] = index[j-inc];
335 }
else if (swapFirst && !swapSecond) {
336 for (
Int_t k = inc; k <
n; k++ ) {
337 const Int_t ktemp = index[k];
341 for (j = k; j >= inc; j -= inc) {
342 if ( fi <
first[j-inc] || (fi ==
first[j-inc] && se <
second[ index[j-inc]])) {
343 index[j] = index[j-inc];
353 for (
Int_t k = inc; k <
n; k++ ) {
354 const Int_t ktemp = index[k];
358 for (j = k; j >= inc; j -= inc) {
360 index [j] = index [j-inc];
383template<
class Element>
391 Element *elem = GetMatrixArray();
393 for (
Int_t irow = 0; irow < fNrows; irow++) {
394 const Int_t off1 = irow*fNcols;
396 for (
Int_t icol = 0; icol < fNcols; icol++) {
397 elem[off1+icol] =
data[off2+irow];
403 memcpy(elem,
data,fNelems*
sizeof(Element));
411template<
class Element>
416 if ((fNrows != fNcols) || (fRowLwb != fColLwb))
419 const Element *
const elem = GetMatrixArray();
420 for (
Int_t irow = 0; irow < fNrows; irow++) {
421 const Int_t rowOff = irow*fNcols;
423 for (
Int_t icol = 0; icol < irow; icol++) {
424 if (elem[rowOff+icol] != elem[colOff+irow])
440template<
class Element>
448 const Element *
const elem = GetMatrixArray();
450 for (
Int_t irow = 0; irow < fNrows; irow++) {
451 const Int_t off1 = irow*fNcols;
453 for (
Int_t icol = 0; icol < fNcols; icol++) {
454 data[off2+irow] = elem[off1+icol];
460 memcpy(
data,elem,fNelems*
sizeof(Element));
466template<
class Element>
469 const Int_t arown = rown-fRowLwb;
470 const Int_t acoln = coln-fColLwb;
471 const Int_t nr = (
n > 0) ?
n : fNcols;
474 if (arown >= fNrows || arown < 0) {
475 Error(
"InsertRow",
"row %d out of matrix range",rown);
479 if (acoln >= fNcols || acoln < 0) {
480 Error(
"InsertRow",
"column %d out of matrix range",coln);
484 if (acoln+nr > fNcols || nr < 0) {
485 Error(
"InsertRow",
"row length %d out of range",nr);
490 const Int_t off = arown*fNcols+acoln;
491 Element *
const elem = GetMatrixArray()+off;
492 memcpy(elem,
v,nr*
sizeof(Element));
500template<
class Element>
503 const Int_t arown = rown-fRowLwb;
504 const Int_t acoln = coln-fColLwb;
505 const Int_t nr = (
n > 0) ?
n : fNcols;
508 if (arown >= fNrows || arown < 0) {
509 Error(
"ExtractRow",
"row %d out of matrix range",rown);
513 if (acoln >= fNcols || acoln < 0) {
514 Error(
"ExtractRow",
"column %d out of matrix range",coln);
518 if (acoln+
n >= fNcols || nr < 0) {
519 Error(
"ExtractRow",
"row length %d out of range",nr);
524 const Int_t off = arown*fNcols+acoln;
525 const Element *
const elem = GetMatrixArray()+off;
526 memcpy(
v,elem,nr*
sizeof(Element));
534template<
class Element>
537 fRowLwb += row_shift;
538 fColLwb += col_shift;
546template<
class Element>
550 memset(this->GetMatrixArray(),0,fNelems*
sizeof(Element));
558template<
class Element>
563 Element *ep = this->GetMatrixArray();
564 const Element *
const fp = ep+fNelems;
576template<
class Element>
581 Element *ep = this->GetMatrixArray();
582 const Element *
const fp = ep+fNelems;
594template<
class Element>
599 Element *ep = this->GetMatrixArray();
600 const Element *
const fp = ep+fNelems;
612template<
class Element>
617 Element *ep = this->GetMatrixArray();
618 memset(ep,0,fNelems*
sizeof(Element));
619 for (
Int_t i = fRowLwb; i <= fRowLwb+fNrows-1; i++)
620 for (
Int_t j = fColLwb; j <= fColLwb+fNcols-1; j++)
621 *ep++ = (i==j ? 1.0 : 0.0);
631template<
class Element>
639 if (
v.GetNoElements() < nMax) {
640 Error(
"NormByDiag",
"vector shorter than matrix diagonal");
649 const Element *pV =
v.GetMatrixArray();
650 Element *mp = this->GetMatrixArray();
653 for (
Int_t irow = 0; irow < fNrows; irow++) {
654 if (pV[irow] != 0.0) {
655 for (
Int_t icol = 0; icol < fNcols; icol++) {
656 if (pV[icol] != 0.0) {
660 Error(
"NormbyDiag",
"vector element %d is zero",icol);
665 Error(
"NormbyDiag",
"vector element %d is zero",irow);
670 for (
Int_t irow = 0; irow < fNrows; irow++) {
671 for (
Int_t icol = 0; icol < fNcols; icol++) {
685template<
class Element>
690 const Element * ep = GetMatrixArray();
691 const Element *
const fp = ep+fNelems;
698 for (
Int_t j = 0; j < fNcols; j++)
712template<
class Element>
717 const Element * ep = GetMatrixArray();
718 const Element *
const fp = ep+fNcols;
725 for (
Int_t i = 0; i < fNrows; i++,ep += fNcols)
739template<
class Element>
744 const Element * ep = GetMatrixArray();
745 const Element *
const fp = ep+fNelems;
748 for ( ; ep < fp; ep++)
749 sum += (*ep) * (*ep);
757template<
class Element>
762 Int_t nr_nonzeros = 0;
763 const Element *ep = this->GetMatrixArray();
764 const Element *
const fp = ep+fNelems;
766 if (*ep++ != 0.0) nr_nonzeros++;
774template<
class Element>
780 const Element *ep = this->GetMatrixArray();
781 const Element *
const fp = ep+fNelems;
791template<
class Element>
796 const Element *
const ep = this->GetMatrixArray();
804template<
class Element>
809 const Element *
const ep = this->GetMatrixArray();
818template<
class Element>
821 gROOT->ProcessLine(
Form(
"THistPainter::PaintSpecialObjects((TObject*)0x%lx,\"%s\");",
831template<
class Element>
835 Error(
"Print",
"Matrix is invalid");
840 const char *format =
"%11.4g ";
842 const char *
f = strstr(option,
"f=");
846 snprintf(topbar,100,format,123.456789);
847 Int_t nch = strlen(topbar)+1;
848 if (nch > 18) nch = 18;
850 for (
Int_t i = 0; i < nch; i++) ftopbar[i] =
' ';
852 snprintf(ftopbar+nch/2,20-nch/2,
"%s%dd",
"%",nk);
853 Int_t nch2 = strlen(ftopbar);
854 for (
Int_t i = nch2; i < nch; i++) ftopbar[i] =
' ';
858 printf(
"\n%dx%d matrix is as follows",fNrows,fNcols);
860 Int_t cols_per_sheet = 5;
861 if (nch <= 8) cols_per_sheet =10;
862 const Int_t ncols = fNcols;
863 const Int_t nrows = fNrows;
864 const Int_t collwb = fColLwb;
865 const Int_t rowlwb = fRowLwb;
867 for (
Int_t i = 0; i < nk; i++) topbar[i] =
'-';
869 for (
Int_t sheet_counter = 1; sheet_counter <= ncols; sheet_counter += cols_per_sheet) {
871 for (
Int_t j = sheet_counter; j < sheet_counter+cols_per_sheet && j <= ncols; j++)
872 printf(ftopbar,j+collwb-1);
873 printf(
"\n%s\n",topbar);
874 if (fNelems <= 0)
continue;
875 for (
Int_t i = 1; i <= nrows; i++) {
876 printf(
"%4d |",i+rowlwb-1);
877 for (
Int_t j = sheet_counter; j < sheet_counter+cols_per_sheet && j <= ncols; j++)
878 printf(format,(*
this)(i+rowlwb-1,j+collwb-1));
888template<
class Element>
893 if (val == 0. && fNelems == 0)
896 const Element * ep = GetMatrixArray();
897 const Element *
const fp = ep+fNelems;
898 for (; ep < fp; ep++)
908template<
class Element>
913 if (val == 0. && fNelems == 0)
916 const Element * ep = GetMatrixArray();
917 const Element *
const fp = ep+fNelems;
918 for (; ep < fp; ep++)
928template<
class Element>
933 const Element * ep = GetMatrixArray();
934 const Element *
const fp = ep+fNelems;
935 for (; ep < fp; ep++)
945template<
class Element>
950 const Element * ep = GetMatrixArray();
951 const Element *
const fp = ep+fNelems;
952 for (; ep < fp; ep++)
962template<
class Element>
967 const Element * ep = GetMatrixArray();
968 const Element *
const fp = ep+fNelems;
969 for (; ep < fp; ep++)
979template<
class Element>
984 const Element * ep = GetMatrixArray();
985 const Element *
const fp = ep+fNelems;
986 for (; ep < fp; ep++)
996template<
class Element>
1001 Element *ep = this->GetMatrixArray();
1002 const Element *
const ep_last = ep+fNelems;
1003 while (ep < ep_last)
1013template<
class Element>
1018 Element *ep = this->GetMatrixArray();
1019 for (action.
fI = fRowLwb; action.
fI < fRowLwb+fNrows; action.
fI++)
1020 for (action.
fJ = fColLwb; action.
fJ < fColLwb+fNcols; action.
fJ++)
1023 R__ASSERT(ep == this->GetMatrixArray()+fNelems);
1031template<
class Element>
1036 const Element scale =
beta-alpha;
1037 const Element shift = alpha/scale;
1039 Element * ep = GetMatrixArray();
1040 const Element *
const fp = ep+fNelems;
1042 *ep++ = scale*(
Drand(seed)+shift);
1050template<
class Element>
1061template<
class Element>
1065 ::Error(
"E2Norm",
"matrices not compatible");
1070 const Element * mp2 =
m2.GetMatrixArray();
1074 for (; mp1 < fmp1; mp1++, mp2++)
1075 sum += (*mp1 - *mp2)*(*mp1 - *mp2);
1083template<
class Element1,
class Element2>
1088 ::Error(
"AreCompatible",
"matrix 1 not valid");
1091 if (!
m2.IsValid()) {
1093 ::Error(
"AreCompatible",
"matrix 2 not valid");
1100 ::Error(
"AreCompatible",
"matrices 1 and 2 not compatible");
1110template<
class Element>
1114 Error(
"Compare(const TMatrixTBase<Element> &,const TMatrixTBase<Element> &)",
"matrices are incompatible");
1118 printf(
"\n\nComparison of two TMatrices:\n");
1125 Element difmax = -1;
1129 const Element mv1 = m1(i,j);
1130 const Element mv2 =
m2(i,j);
1133 if (diff > difmax) {
1144 printf(
"\nMaximal discrepancy \t\t%g", difmax);
1145 printf(
"\n occured at the point\t\t(%d,%d)",imax,jmax);
1146 const Element mv1 = m1(imax,jmax);
1147 const Element mv2 =
m2(imax,jmax);
1148 printf(
"\n Matrix 1 element is \t\t%g", mv1);
1149 printf(
"\n Matrix 2 element is \t\t%g", mv2);
1150 printf(
"\n Absolute error v2[i]-v1[i]\t\t%g", mv2-mv1);
1151 printf(
"\n Relative error\t\t\t\t%g\n",
1154 printf(
"\n||Matrix 1|| \t\t\t%g", norm1);
1155 printf(
"\n||Matrix 2|| \t\t\t%g", norm2);
1156 printf(
"\n||Matrix1-Matrix2||\t\t\t\t%g", ndiff);
1157 printf(
"\n||Matrix1-Matrix2||/sqrt(||Matrix1|| ||Matrix2||)\t%g\n\n",
1164template<
class Element>
1174 Element maxDevObs = 0;
1179 for (
Int_t i =
m.GetRowLwb(); i <=
m.GetRowUpb(); i++) {
1180 for (
Int_t j =
m.GetColLwb(); j <=
m.GetColUpb(); j++) {
1182 if (dev > maxDevObs) {
1194 printf(
"Largest dev for (%d,%d); dev = |%g - %g| = %g\n",imax,jmax,
m(imax,jmax),val,maxDevObs);
1195 if(maxDevObs > maxDevAllow)
1196 Error(
"VerifyElementValue",
"Deviation > %g\n",maxDevAllow);
1199 if(maxDevObs > maxDevAllow)
1207template<
class Element>
1209 Element maxDevAllow)
1214 if (m1 == 0 &&
m2 == 0)
1219 Element maxDevObs = 0;
1227 if (dev > maxDevObs) {
1239 printf(
"Largest dev for (%d,%d); dev = |%g - %g| = %g\n",
1240 imax,jmax,m1(imax,jmax),
m2(imax,jmax),maxDevObs);
1241 if (maxDevObs > maxDevAllow)
1242 Error(
"VerifyMatrixValue",
"Deviation > %g\n",maxDevAllow);
1245 if (maxDevObs > maxDevAllow)
1253template<
class Element>
1262 Error(
"TMatrixTBase<Element>::Streamer",
"Unknown version number: %d",R__v);
1265 if (R__v < 4) MakeValid();
1272template<
class Element>
1274 static Element gNanValue;
1277Double_t nan_value_t<Double_t>::gNanValue = std::numeric_limits<Double_t>::quiet_NaN();
1279Float_t nan_value_t<Float_t>::gNanValue = std::numeric_limits<Float_t>::quiet_NaN();
1281template<
class Element>
1284 return nan_value_t<Element>::gNanValue;
#define templateClassImp(name)
void Error(const char *location, const char *msgfmt,...)
template void Compare< Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2)
template Bool_t VerifyMatrixValue< Double_t >(const TMatrixDBase &m, Double_t val, Int_t verbose, Double_t maxDevAllow)
template void Compare< Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2)
Bool_t VerifyMatrixIdentity(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2, Int_t verbose, Element maxDevAllow)
Verify that elements of the two matrices are equal within MaxDevAllow .
template Bool_t AreCompatible< Float_t, Double_t >(const TMatrixFBase &m1, const TMatrixDBase &m2, Int_t verbose)
template Bool_t AreCompatible< Float_t, Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2, Int_t verbose)
void Compare(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Compare two matrices and print out the result of the comparison.
Element E2Norm(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Square of the Euclidian norm of the difference between two matrices.
template Bool_t VerifyMatrixIdentity< Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2, Int_t verbose, Double_t maxDevAllow)
Int_t gMatrixCheck
TMatrixTBase.
template Double_t E2Norm< Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2)
Bool_t VerifyMatrixValue(const TMatrixTBase< Element > &m, Element val, Int_t verbose, Element maxDevAllow)
Validate that all elements of matrix have value val within maxDevAllow.
template Float_t E2Norm< Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2)
template Bool_t AreCompatible< Double_t, Double_t >(const TMatrixDBase &m1, const TMatrixDBase &m2, Int_t verbose)
template Bool_t VerifyMatrixIdentity< Float_t >(const TMatrixFBase &m1, const TMatrixFBase &m2, Int_t verbose, Float_t maxDevAllowN)
template Bool_t AreCompatible< Double_t, Float_t >(const TMatrixDBase &m1, const TMatrixFBase &m2, Int_t verbose)
Bool_t operator==(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Check to see if two matrices are identical.
template Bool_t VerifyMatrixValue< Float_t >(const TMatrixFBase &m, Float_t val, Int_t verbose, Float_t maxDevAllow)
Bool_t AreCompatible(const TMatrixTBase< Element1 > &m1, const TMatrixTBase< Element2 > &m2, Int_t verbose)
Check that matrice sm1 and m2 areboth valid and have identical shapes .
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
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
double beta(double x, double y)
Calculates the beta function.
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.
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)
static long int sum(long int i)