46 template<
class Element>
49 const Int_t no_rows = row_upb-row_lwb+1;
50 Allocate(no_rows,no_rows,row_lwb,row_lwb,1);
61 template<
class Element>
65 SetMatrixArray(elements,option);
66 if (!this->IsSymmetric()) {
67 Error(
"TMatrixTSym(Int_t,Element*,Option_t*)",
"matrix not symmetric");
74 template<
class Element>
77 const Int_t no_rows = row_upb-row_lwb+1;
78 Allocate(no_rows,no_rows,row_lwb,row_lwb);
79 SetMatrixArray(elements,option);
80 if (!this->IsSymmetric()) {
81 Error(
"TMatrixTSym(Int_t,Int_t,Element*,Option_t*)",
"matrix not symmetric");
87 template<
class Element>
100 template<
class Element>
132 this->SetTol(oldTol);
142 Error(
"TMatrixTSym(EMatrixCreatorOp1,const TMatrixTSym)",
143 "operation %d not yet implemented", op);
149 template<
class Element>
161 Error(
"TMatrixTSym(EMatrixCreatorOp1,const TMatrixT)",
162 "operation %d not yet implemented", op);
168 template<
class Element>
190 Error(
"TMatrixTSym(EMatrixCreatorOp2)",
"operation %d not yet implemented", op);
196 template<
class Element>
202 lazy_constructor.
FillIn(*
this);
203 if (!this->IsSymmetric()) {
204 Error(
"TMatrixTSym(TMatrixTSymLazy)",
"matrix not symmetric");
211 template<
class Element>
215 if (size > this->kSizeMax)
225 template<
class Element>
228 if (size == 0)
return 0;
230 if ( size <= this->kSizeMax )
233 Element *heap =
new Element[size];
243 template<
class Element>
247 if (copySize == 0 || oldp == newp)
250 if ( newSize <= this->kSizeMax && oldSize <= this->kSizeMax ) {
253 for (
Int_t i = copySize-1; i >= 0; i--)
256 for (
Int_t i = 0; i < copySize; i++)
261 memcpy(newp,oldp,copySize*
sizeof(Element));
270 template<
class Element>
274 this->fIsOwner =
kTRUE;
283 if (no_rows < 0 || no_cols < 0)
285 Error(
"Allocate",
"no_rows=%d no_cols=%d",no_rows,no_cols);
291 this->fNrows = no_rows;
292 this->fNcols = no_cols;
293 this->fRowLwb = row_lwb;
294 this->fColLwb = col_lwb;
295 this->fNelems = this->fNrows*this->fNcols;
297 if (this->fNelems > 0) {
298 fElements = New_m(this->fNelems);
300 memset(fElements,0,this->fNelems*
sizeof(Element));
308 template<
class Element>
313 Error(
"Plus",
"matrices not compatible");
318 Error(
"Plus",
"this->GetMatrixArray() == a.GetMatrixArray()");
323 Error(
"Plus",
"this->GetMatrixArray() == b.GetMatrixArray()");
330 Element * cp = this->GetMatrixArray();
331 const Element *
const cp_last = cp+this->fNelems;
333 while (cp < cp_last) {
342 template<
class Element>
347 Error(
"Minus",
"matrices not compatible");
352 Error(
"Minus",
"this->GetMatrixArray() == a.GetMatrixArray()");
357 Error(
"Minus",
"this->GetMatrixArray() == b.GetMatrixArray()");
364 Element * cp = this->GetMatrixArray();
365 const Element *
const cp_last = cp+this->fNelems;
367 while (cp < cp_last) {
377 template<
class Element>
384 Element *cp = this->GetMatrixArray();
385 if (
typeid(Element) ==
typeid(
Double_t))
386 cblas_dgemm (CblasRowMajor,CblasTrans,CblasNoTrans,this->fNrows,this->fNcols,a.
GetNrows(),
388 else if (
typeid(Element) !=
typeid(
Float_t))
389 cblas_sgemm (CblasRowMajor,CblasTrans,CblasNoTrans,fNrows,fNcols,a.
GetNrows(),
392 Error(
"TMult",
"type %s not implemented in BLAS library",
typeid(Element));
396 const Int_t ncolsb = ncolsa;
398 const Element *
const bp = ap;
399 Element * cp = this->GetMatrixArray();
401 const Element *acp0 = ap;
403 for (
const Element *bcp = bp; bcp < bp+ncolsb; ) {
404 const Element *acp = acp0;
406 while (bcp < bp+nb) {
417 R__ASSERT(cp == this->GetMatrixArray()+this->fNelems && acp0 == ap+ncolsa);
425 template<
class Element>
432 Element *cp = this->GetMatrixArray();
433 if (
typeid(Element) ==
typeid(
Double_t))
434 cblas_dsymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
436 else if (
typeid(Element) !=
typeid(
Float_t))
437 cblas_ssymm (CblasRowMajor,CblasLeft,CblasUpper,fNrows,fNcols,1.0,
440 Error(
"TMult",
"type %s not implemented in BLAS library",
typeid(Element));
444 const Int_t ncolsb = ncolsa;
446 const Element *
const bp = ap;
447 Element * cp = this->GetMatrixArray();
449 const Element *acp0 = ap;
451 for (
const Element *bcp = bp; bcp < bp+ncolsb; ) {
452 const Element *acp = acp0;
454 while (bcp < bp+nb) {
465 R__ASSERT(cp == this->GetMatrixArray()+this->fNelems && acp0 == ap+ncolsa);
471 template<
class Element>
476 Error(
"Use",
"row_upb=%d < row_lwb=%d",row_upb,row_lwb);
481 this->fNrows = row_upb-row_lwb+1;
482 this->fNcols = this->fNrows;
483 this->fRowLwb = row_lwb;
484 this->fColLwb = row_lwb;
485 this->fNelems = this->fNrows*this->fNcols;
499 template<
class Element>
505 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
506 Error(
"GetSub",
"row_lwb out of bounds");
509 if (row_upb < this->fRowLwb || row_upb > this->fRowLwb+this->fNrows-1) {
510 Error(
"GetSub",
"row_upb out of bounds");
513 if (row_upb < row_lwb) {
514 Error(
"GetSub",
"row_upb < row_lwb");
527 row_upb_sub = row_upb-row_lwb;
529 row_lwb_sub = row_lwb;
530 row_upb_sub = row_upb;
533 target.
ResizeTo(row_lwb_sub,row_upb_sub,row_lwb_sub,row_upb_sub);
534 const Int_t nrows_sub = row_upb_sub-row_lwb_sub+1;
537 for (
Int_t irow = 0; irow < nrows_sub; irow++) {
538 for (
Int_t icol = 0; icol < nrows_sub; icol++) {
539 target(irow+row_lwb_sub,icol+row_lwb_sub) = (*this)(row_lwb+irow,row_lwb+icol);
543 const Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNrows+(row_lwb-this->fRowLwb);
546 for (
Int_t irow = 0; irow < nrows_sub; irow++) {
547 const Element *ap_sub = ap;
548 for (
Int_t icol = 0; icol < nrows_sub; icol++) {
565 template<
class Element>
571 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
572 Error(
"GetSub",
"row_lwb out of bounds");
575 if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
576 Error(
"GetSub",
"col_lwb out of bounds");
579 if (row_upb < this->fRowLwb || row_upb > this->fRowLwb+this->fNrows-1) {
580 Error(
"GetSub",
"row_upb out of bounds");
583 if (col_upb < this->fColLwb || col_upb > this->fColLwb+this->fNcols-1) {
584 Error(
"GetSub",
"col_upb out of bounds");
587 if (row_upb < row_lwb || col_upb < col_lwb) {
588 Error(
"GetSub",
"row_upb < row_lwb || col_upb < col_lwb");
597 const Int_t row_lwb_sub = (shift) ? 0 : row_lwb;
598 const Int_t row_upb_sub = (shift) ? row_upb-row_lwb : row_upb;
599 const Int_t col_lwb_sub = (shift) ? 0 : col_lwb;
600 const Int_t col_upb_sub = (shift) ? col_upb-col_lwb : col_upb;
602 target.
ResizeTo(row_lwb_sub,row_upb_sub,col_lwb_sub,col_upb_sub);
603 const Int_t nrows_sub = row_upb_sub-row_lwb_sub+1;
604 const Int_t ncols_sub = col_upb_sub-col_lwb_sub+1;
607 for (
Int_t irow = 0; irow < nrows_sub; irow++) {
608 for (
Int_t icol = 0; icol < ncols_sub; icol++) {
609 target(irow+row_lwb_sub,icol+col_lwb_sub) = (*this)(row_lwb+irow,col_lwb+icol);
613 const Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNcols+(col_lwb-this->fColLwb);
616 for (
Int_t irow = 0; irow < nrows_sub; irow++) {
617 const Element *ap_sub = ap;
618 for (
Int_t icol = 0; icol < ncols_sub; icol++) {
632 template<
class Element>
640 Error(
"SetSub",
"source matrix is not symmetric");
643 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
644 Error(
"SetSub",
"row_lwb outof bounds");
647 if (row_lwb+source.
GetNrows() > this->fRowLwb+this->fNrows) {
648 Error(
"SetSub",
"source matrix too large");
657 for (
Int_t irow = 0; irow < nRows_source; irow++) {
658 for (
Int_t icol = 0; icol < nRows_source; icol++) {
659 (*this)(row_lwb+irow,row_lwb+icol) = source(rowlwb_s+irow,rowlwb_s+icol);
664 Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNrows+(row_lwb-this->fRowLwb);
666 for (
Int_t irow = 0; irow < nRows_source; irow++) {
667 Element *ap_sub = ap;
668 for (
Int_t icol = 0; icol < nRows_source; icol++) {
682 template<
class Element>
689 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
690 Error(
"SetSub",
"row_lwb out of bounds");
693 if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
694 Error(
"SetSub",
"col_lwb out of bounds");
698 if (row_lwb+source.
GetNrows() > this->fRowLwb+this->fNrows || col_lwb+source.
GetNcols() > this->fRowLwb+this->fNrows) {
699 Error(
"SetSub",
"source matrix too large");
702 if (col_lwb+source.
GetNcols() > this->fRowLwb+this->fNrows || row_lwb+source.
GetNrows() > this->fRowLwb+this->fNrows) {
703 Error(
"SetSub",
"source matrix too large");
713 if (row_lwb >= col_lwb) {
716 for (irow = 0; irow < nRows_source; irow++) {
717 for (
Int_t icol = 0; col_lwb+icol <= row_lwb+irow &&
718 icol < nCols_source; icol++) {
719 (*this)(row_lwb+irow,col_lwb+icol) = source(irow+rowlwb_s,icol+collwb_s);
724 for (irow = 0; irow < nCols_source; irow++) {
725 for (
Int_t icol = nRows_source-1; row_lwb+icol > irow+col_lwb &&
727 (*this)(col_lwb+irow,row_lwb+icol) = source(icol+rowlwb_s,irow+collwb_s);
739 template<
class Element>
743 if (!this->IsSymmetric()) {
744 Error(
"SetMatrixArray",
"Matrix is not symmetric after Set");
752 template<
class Element>
755 if (row_shift != col_shift) {
756 Error(
"Shift",
"row_shift != col_shift");
767 template<
class Element>
771 if (!this->fIsOwner) {
772 Error(
"ResizeTo(Int_t,Int_t)",
"Not owner of data array,cannot resize");
776 if (nrows != ncols) {
777 Error(
"ResizeTo(Int_t,Int_t)",
"nrows != ncols");
781 if (this->fNelems > 0) {
782 if (this->fNrows == nrows && this->fNcols == ncols)
784 else if (nrows == 0 || ncols == 0) {
785 this->fNrows = nrows; this->fNcols = ncols;
790 Element *elements_old = GetMatrixArray();
791 const Int_t nelems_old = this->fNelems;
792 const Int_t nrows_old = this->fNrows;
793 const Int_t ncols_old = this->fNcols;
798 Element *elements_new = GetMatrixArray();
801 if (this->fNelems > this->kSizeMax || nelems_old > this->kSizeMax)
802 memset(elements_new,0,this->fNelems*
sizeof(Element));
803 else if (this->fNelems > nelems_old)
804 memset(elements_new+nelems_old,0,(this->fNelems-nelems_old)*
sizeof(Element));
810 const Int_t nelems_new = this->fNelems;
811 if (ncols_old < this->fNcols) {
812 for (
Int_t i = nrows_copy-1; i >= 0; i--) {
813 Memcpy_m(elements_new+i*this->fNcols,elements_old+i*ncols_old,ncols_copy,
814 nelems_new,nelems_old);
815 if (this->fNelems <= this->kSizeMax && nelems_old <= this->kSizeMax)
816 memset(elements_new+i*this->fNcols+ncols_copy,0,(this->fNcols-ncols_copy)*
sizeof(Element));
819 for (
Int_t i = 0; i < nrows_copy; i++)
820 Memcpy_m(elements_new+i*this->fNcols,elements_old+i*ncols_old,ncols_copy,
821 nelems_new,nelems_old);
824 Delete_m(nelems_old,elements_old);
837 template<
class Element>
842 if (!this->fIsOwner) {
843 Error(
"ResizeTo(Int_t,Int_t,Int_t,Int_t)",
"Not owner of data array,cannot resize");
847 if (row_lwb != col_lwb) {
848 Error(
"ResizeTo(Int_t,Int_t,Int_t,Int_t)",
"row_lwb != col_lwb");
851 if (row_upb != col_upb) {
852 Error(
"ResizeTo(Int_t,Int_t,Int_t,Int_t)",
"row_upb != col_upb");
856 const Int_t new_nrows = row_upb-row_lwb+1;
857 const Int_t new_ncols = col_upb-col_lwb+1;
859 if (this->fNelems > 0) {
861 if (this->fNrows == new_nrows && this->fNcols == new_ncols &&
862 this->fRowLwb == row_lwb && this->fColLwb == col_lwb)
864 else if (new_nrows == 0 || new_ncols == 0) {
865 this->fNrows = new_nrows; this->fNcols = new_ncols;
866 this->fRowLwb = row_lwb; this->fColLwb = col_lwb;
871 Element *elements_old = GetMatrixArray();
872 const Int_t nelems_old = this->fNelems;
873 const Int_t nrows_old = this->fNrows;
874 const Int_t ncols_old = this->fNcols;
875 const Int_t rowLwb_old = this->fRowLwb;
876 const Int_t colLwb_old = this->fColLwb;
878 Allocate(new_nrows,new_ncols,row_lwb,col_lwb);
881 Element *elements_new = GetMatrixArray();
884 if (this->fNelems > this->kSizeMax || nelems_old > this->kSizeMax)
885 memset(elements_new,0,this->fNelems*
sizeof(Element));
886 else if (this->fNelems > nelems_old)
887 memset(elements_new+nelems_old,0,(this->fNelems-nelems_old)*
sizeof(Element));
892 const Int_t rowUpb_copy =
TMath::Min(this->fRowLwb+this->fNrows-1,rowLwb_old+nrows_old-1);
893 const Int_t colUpb_copy =
TMath::Min(this->fColLwb+this->fNcols-1,colLwb_old+ncols_old-1);
895 const Int_t nrows_copy = rowUpb_copy-rowLwb_copy+1;
896 const Int_t ncols_copy = colUpb_copy-colLwb_copy+1;
898 if (nrows_copy > 0 && ncols_copy > 0) {
899 const Int_t colOldOff = colLwb_copy-colLwb_old;
900 const Int_t colNewOff = colLwb_copy-this->fColLwb;
901 if (ncols_old < this->fNcols) {
902 for (
Int_t i = nrows_copy-1; i >= 0; i--) {
903 const Int_t iRowOld = rowLwb_copy+i-rowLwb_old;
904 const Int_t iRowNew = rowLwb_copy+i-this->fRowLwb;
905 Memcpy_m(elements_new+iRowNew*this->fNcols+colNewOff,
906 elements_old+iRowOld*ncols_old+colOldOff,ncols_copy,this->fNelems,nelems_old);
907 if (this->fNelems <= this->kSizeMax && nelems_old <= this->kSizeMax)
908 memset(elements_new+iRowNew*this->fNcols+colNewOff+ncols_copy,0,
909 (this->fNcols-ncols_copy)*
sizeof(Element));
912 for (
Int_t i = 0; i < nrows_copy; i++) {
913 const Int_t iRowOld = rowLwb_copy+i-rowLwb_old;
914 const Int_t iRowNew = rowLwb_copy+i-this->fRowLwb;
915 Memcpy_m(elements_new+iRowNew*this->fNcols+colNewOff,
916 elements_old+iRowOld*ncols_old+colOldOff,ncols_copy,this->fNelems,nelems_old);
921 Delete_m(nelems_old,elements_old);
923 Allocate(new_nrows,new_ncols,row_lwb,col_lwb,1);
931 template<
class Element>
943 template<
class Element>
957 template<
class Element>
964 Element *
p2 = this->GetMatrixArray();
965 for (
Int_t i = 0; i < this->GetNoElements(); i++)
975 template<
class Element>
984 Element *pM = this->GetMatrixArray();
986 Error(
"InvertFast",
"matrix is singular");
996 TMatrixTSymCramerInv::Inv2x2<Element>(*
this,det);
1001 TMatrixTSymCramerInv::Inv3x3<Element>(*
this,det);
1006 TMatrixTSymCramerInv::Inv4x4<Element>(*
this,det);
1011 TMatrixTSymCramerInv::Inv5x5<Element>(*
this,det);
1016 TMatrixTSymCramerInv::Inv6x6<Element>(*
this,det);
1025 Element *
p2 = this->GetMatrixArray();
1026 for (
Int_t i = 0; i < this->GetNoElements(); i++)
1037 template<
class Element>
1046 Error(
"Transpose",
"matrix has wrong shape");
1059 template<
class Element>
1066 Error(
"Rank1Update",
"vector too short");
1072 Element *trp = this->GetMatrixArray();
1074 for (
Int_t i = 0; i < this->fNrows; i++) {
1076 tcp += i*this->fNcols;
1077 const Element tmp = alpha*pv[i];
1078 for (
Int_t j = i; j < this->fNcols; j++) {
1079 if (j > i) *tcp += tmp*pv[j];
1080 *trp++ += tmp*pv[j];
1081 tcp += this->fNcols;
1083 tcp -= this->fNelems-1;
1095 template<
class Element>
1102 Error(
"Similarity(const TMatrixT &)",
"matrices incompatible");
1107 const Int_t ncolsa = this->fNcols;
1114 Element work[kWorkMax];
1116 Element *bap = work;
1117 if (nrowsb*ncolsa > kWorkMax) {
1118 isAllocated =
kTRUE;
1119 bap =
new Element[nrowsb*ncolsa];
1122 AMultB(bp,nb,ncolsb,this->fElements,this->fNelems,this->fNcols,bap);
1124 if (nrowsb != this->fNrows)
1125 this->ResizeTo(nrowsb,nrowsb);
1128 Element *cp = this->GetMatrixArray();
1129 if (
typeid(Element) ==
typeid(
Double_t))
1130 cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasTrans,this->fNrows,this->fNcols,ba.GetNcols(),
1131 1.0,bap,ba.GetNcols(),bp,b.
GetNcols(),1.0,cp,this->fNcols);
1132 else if (
typeid(Element) !=
typeid(
Float_t))
1133 cblas_sgemm (CblasRowMajor,CblasNoTrans,CblasTrans,this->fNrows,this->fNcols,ba.GetNcols(),
1134 1.0,bap,ba.GetNcols(),bp,b.
GetNcols(),1.0,cp,this->fNcols);
1136 Error(
"Similarity",
"type %s not implemented in BLAS library",
typeid(Element));
1138 const Int_t nba = nrowsb*ncolsa;
1139 const Int_t ncolsba = ncolsa;
1140 const Element * bi1p = bp;
1141 Element * cp = this->GetMatrixArray();
1142 Element *
const cp0 = cp;
1145 const Element *barp0 = bap;
1146 while (barp0 < bap+nba) {
1147 const Element *brp0 = bi1p;
1148 while (brp0 < bp+nb) {
1149 const Element *barp = barp0;
1150 const Element *brp = brp0;
1152 while (brp < brp0+ncolsb)
1153 cij += *barp++ * *brp++;
1162 R__ASSERT(cp == cp0+this->fNelems+ishift && barp0 == bap+nba);
1165 for (
Int_t irow = 0; irow < this->fNrows; irow++) {
1166 const Int_t rowOff1 = irow*this->fNrows;
1167 for (
Int_t icol = 0; icol < irow; icol++) {
1168 const Int_t rowOff2 = icol*this->fNrows;
1169 cp[rowOff1+icol] = cp[rowOff2+irow];
1186 template<
class Element>
1193 Error(
"Similarity(const TMatrixTSym &)",
"matrices incompatible");
1200 const Int_t ncolsa = this->GetNcols();
1202 Element work[kWorkMax];
1204 Element *abtp = work;
1205 if (this->fNcols > kWorkMax) {
1206 isAllocated =
kTRUE;
1207 abtp =
new Element[this->fNcols];
1213 Element *cp = this->GetMatrixArray();
1214 if (
typeid(Element) ==
typeid(
Double_t))
1215 cblas_dsymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
1217 else if (
typeid(Element) !=
typeid(
Float_t))
1218 cblas_ssymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
1221 Error(
"Similarity",
"type %s not implemented in BLAS library",
typeid(Element));
1226 const Int_t ncolsa = this->GetNcols();
1233 Element work[kWorkMax];
1235 Element *bap = work;
1236 if (nrowsb*ncolsa > kWorkMax) {
1237 isAllocated =
kTRUE;
1238 bap =
new Element[nrowsb*ncolsa];
1241 AMultB(bp,nb,ncolsb,this->fElements,this->fNelems,this->fNcols,bap);
1243 const Int_t nba = nrowsb*ncolsa;
1244 const Int_t ncolsba = ncolsa;
1245 const Element * bi1p = bp;
1246 Element * cp = this->GetMatrixArray();
1247 Element *
const cp0 = cp;
1250 const Element *barp0 = bap;
1251 while (barp0 < bap+nba) {
1252 const Element *brp0 = bi1p;
1253 while (brp0 < bp+nb) {
1254 const Element *barp = barp0;
1255 const Element *brp = brp0;
1257 while (brp < brp0+ncolsb)
1258 cij += *barp++ * *brp++;
1267 R__ASSERT(cp == cp0+this->fNelems+ishift && barp0 == bap+nba);
1270 for (
Int_t irow = 0; irow < this->fNrows; irow++) {
1271 const Int_t rowOff1 = irow*this->fNrows;
1272 for (
Int_t icol = 0; icol < irow; icol++) {
1273 const Int_t rowOff2 = icol*this->fNrows;
1274 cp[rowOff1+icol] = cp[rowOff2+irow];
1288 template<
class Element>
1294 if (this->fNcols != v.
GetNrows() || this->fColLwb != v.
GetLwb()) {
1295 Error(
"Similarity(const TVectorT &)",
"vector and matrix incompatible");
1300 const Element *mp = this->GetMatrixArray();
1304 const Element *
const vp_first = vp;
1305 const Element *
const vp_last = vp+v.
GetNrows();
1306 while (vp < vp_last) {
1308 for (
const Element *sp = vp_first; sp < vp_last; )
1309 sum2 += *mp++ * *sp++;
1310 sum1 += sum2 * *vp++;
1313 R__ASSERT(mp == this->GetMatrixArray()+this->GetNoElements());
1323 template<
class Element>
1330 Error(
"SimilarityT(const TMatrixT &)",
"matrices incompatible");
1336 const Int_t ncolsa = this->GetNcols();
1338 Element work[kWorkMax];
1340 Element *btap = work;
1341 if (ncolsb*ncolsa > kWorkMax) {
1342 isAllocated =
kTRUE;
1343 btap =
new Element[ncolsb*ncolsa];
1349 if (ncolsb != this->fNcols)
1350 this->ResizeTo(ncolsb,ncolsb);
1354 Element *cp = this->GetMatrixArray();
1355 if (
typeid(Element) ==
typeid(
Double_t))
1356 cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasNoTrans,this->fNrows,this->fNcols,bta.
GetNcols(),
1358 else if (
typeid(Element) !=
typeid(
Float_t))
1359 cblas_sgemm (CblasRowMajor,CblasNoTrans,CblasNoTrans,this->fNrows,this->fNcols,bta.
GetNcols(),
1362 Error(
"similarityT",
"type %s not implemented in BLAS library",
typeid(Element));
1368 Element * cp = this->GetMatrixArray();
1369 Element *
const cp0 = cp;
1372 const Element *btarp0 = btap;
1373 const Element *bcp0 = bp;
1374 while (btarp0 < btap+nbta) {
1375 for (
const Element *bcp = bcp0; bcp < bp+ncolsb; ) {
1376 const Element *btarp = btarp0;
1378 while (bcp < bp+nb) {
1379 cij += *btarp++ * *bcp;
1390 R__ASSERT(cp == cp0+this->fNelems+ishift && btarp0 == btap+nbta);
1393 for (
Int_t irow = 0; irow < this->fNrows; irow++) {
1394 const Int_t rowOff1 = irow*this->fNrows;
1395 for (
Int_t icol = 0; icol < irow; icol++) {
1396 const Int_t rowOff2 = icol*this->fNrows;
1397 cp[rowOff1+icol] = cp[rowOff2+irow];
1410 template<
class Element>
1414 Error(
"operator=",
"matrices not compatible");
1420 memcpy(this->GetMatrixArray(),source.
fElements,this->fNelems*
sizeof(Element));
1427 template<
class Element>
1432 if (lazy_constructor.
fRowUpb != this->GetRowUpb() ||
1433 lazy_constructor.
fRowLwb != this->GetRowLwb()) {
1434 Error(
"operator=(const TMatrixTSymLazy&)",
"matrix is incompatible with "
1435 "the assigned Lazy matrix");
1439 lazy_constructor.
FillIn(*
this);
1446 template<
class Element>
1451 Element *ep = fElements;
1452 const Element *
const ep_last = ep+this->fNelems;
1453 while (ep < ep_last)
1462 template<
class Element>
1467 Element *ep = fElements;
1468 const Element *
const ep_last = ep+this->fNelems;
1469 while (ep < ep_last)
1478 template<
class Element>
1483 Element *ep = fElements;
1484 const Element *
const ep_last = ep+this->fNelems;
1485 while (ep < ep_last)
1494 template<
class Element>
1499 Element *ep = fElements;
1500 const Element *
const ep_last = ep+this->fNelems;
1501 while (ep < ep_last)
1510 template<
class Element>
1514 Error(
"operator+=",
"matrices not compatible");
1519 Element *tp = this->GetMatrixArray();
1520 const Element *
const tp_last = tp+this->fNelems;
1521 while (tp < tp_last)
1530 template<
class Element>
1534 Error(
"operator-=",
"matrices not compatible");
1539 Element *tp = this->GetMatrixArray();
1540 const Element *
const tp_last = tp+this->fNelems;
1541 while (tp < tp_last)
1549 template<
class Element>
1555 Element *trp = this->GetMatrixArray();
1557 for (
Int_t i = 0; i < this->fNrows; i++) {
1559 tcp += i*this->fNcols;
1560 for (
Int_t j = i; j < this->fNcols; j++) {
1562 if (j > i) *tcp = val;
1564 tcp += this->fNcols;
1566 tcp -= this->fNelems-1;
1576 template<
class Element>
1582 Element *trp = this->GetMatrixArray();
1584 for (
Int_t i = 0; i < this->fNrows; i++) {
1585 action.
fI = i+this->fRowLwb;
1587 tcp += i*this->fNcols;
1588 for (
Int_t j = i; j < this->fNcols; j++) {
1589 action.
fJ = j+this->fColLwb;
1591 if (j > i) *tcp = val;
1593 tcp += this->fNcols;
1595 tcp -= this->fNelems-1;
1604 template<
class Element>
1609 if (this->fNrows != this->fNcols || this->fRowLwb != this->fColLwb) {
1610 Error(
"Randomize(Element,Element,Element &",
"matrix should be square");
1615 const Element scale = beta-alpha;
1616 const Element shift = alpha/scale;
1618 Element *ep = GetMatrixArray();
1619 for (
Int_t i = 0; i < this->fNrows; i++) {
1620 const Int_t off = i*this->fNcols;
1621 for (
Int_t j = 0; j <= i; j++) {
1622 ep[off+j] = scale*(
Drand(seed)+shift);
1624 ep[j*this->fNcols+i] = ep[off+j];
1635 template<
class Element>
1640 if (this->fNrows != this->fNcols || this->fRowLwb != this->fColLwb) {
1641 Error(
"RandomizeSym(Element,Element,Element &",
"matrix should be square");
1646 const Element scale = beta-alpha;
1647 const Element shift = alpha/scale;
1649 Element *ep = GetMatrixArray();
1651 for (i = 0; i < this->fNrows; i++) {
1652 const Int_t off = i*this->fNcols;
1653 for (
Int_t j = 0; j <= i; j++)
1654 ep[off+j] = scale*(
Drand(seed)+shift);
1657 for (i = this->fNrows-1; i >= 0; i--) {
1658 const Int_t off1 = i*this->fNcols;
1659 for (
Int_t j = i; j >= 0; j--) {
1660 const Int_t off2 = j*this->fNcols;
1661 ep[off1+j] *= ep[off2+j];
1662 for (
Int_t k = j-1; k >= 0; k--) {
1663 ep[off1+j] += ep[off1+k]*ep[off2+k];
1666 ep[off2+i] = ep[off1+j];
1677 template<
class Element>
1680 TMatrixDSym tmp = *
this;
1682 eigenValues.
ResizeTo(this->fNrows);
1690 template<
class Element>
1700 template<
class Element>
1710 template<
class Element>
1720 template<
class Element>
1728 template<
class Element>
1738 template<
class Element>
1748 template<
class Element>
1751 return Element(-1.0)*
operator-(source1,val);
1756 template<
class Element>
1766 template<
class Element>
1775 template<
class Element>
1781 Error(
"operator&&(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1791 while (tp < tp_last)
1792 *tp++ = (*sp1++ != 0.0 && *sp2++ != 0.0);
1800 template<
class Element>
1806 Error(
"operator||(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1816 while (tp < tp_last)
1817 *tp++ = (*sp1++ != 0.0 || *sp2++ != 0.0);
1825 template<
class Element>
1831 Error(
"operator>(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1841 while (tp < tp_last) {
1842 *tp++ = (*sp1) > (*sp2); sp1++; sp2++;
1851 template<
class Element>
1857 Error(
"operator>=(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1867 while (tp < tp_last) {
1868 *tp++ = (*sp1) >= (*sp2); sp1++; sp2++;
1877 template<
class Element>
1883 Error(
"operator<=(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1890 const Element *sp2 = source2.GetMatrixArray();
1893 while (tp < tp_last) {
1894 *tp++ = (*sp1) <= (*sp2); sp1++; sp2++;
1903 template<
class Element>
1909 Error(
"operator<(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1916 const Element *sp2 = source2.GetMatrixArray();
1919 while (tp < tp_last) {
1920 *tp++ = (*sp1) < (*sp2); sp1++; sp2++;
1929 template<
class Element>
1933 ::Error(
"Add",
"matrices not compatible");
1943 for (
Int_t i = 0; i < nrows; i++) {
1947 for (
Int_t j = i; j < ncols; j++) {
1948 const Element tmp = scalar * *sp++;
1949 if (j > i) *tcp += tmp;
1962 template<
class Element>
1966 ::Error(
"ElementMult",
"matrices not compatible");
1976 for (
Int_t i = 0; i < nrows; i++) {
1980 for (
Int_t j = i; j < ncols; j++) {
1981 if (j > i) *tcp *= *sp;
1994 template<
class Element>
1998 ::Error(
"ElementDiv",
"matrices not compatible");
2008 for (
Int_t i = 0; i < nrows; i++) {
2012 for (
Int_t j = i; j < ncols; j++) {
2014 if (j > i) *tcp /= *sp;
2019 Error(
"ElementDiv",
"source (%d,%d) is zero",irow,icol);
2033 template<
class Element>
2041 fElements =
new Element[this->fNelems];
2043 for (i = 0; i < this->fNrows; i++) {
2044 R__b.
ReadFastArray(fElements+i*this->fNcols+i,this->fNcols-i);
2047 for (i = 0; i < this->fNrows; i++) {
2048 for (
Int_t j = 0; j < i; j++) {
2049 fElements[i*this->fNcols+j] = fElements[j*this->fNrows+i];
2052 if (this->fNelems <= this->kSizeMax) {
2053 memcpy(fDataStack,fElements,this->fNelems*
sizeof(Element));
2054 delete [] fElements;
2055 fElements = fDataStack;
2060 for (
Int_t i = 0; i < this->fNrows; i++) {
2066 #ifndef ROOT_TMatrixFSymfwd
2072 template Bool_t operator== <
Float_t>(
const TMatrixFSym &source1,
const TMatrixFSym &source2);
2073 template TMatrixFSym
operator+ <
Float_t>(
const TMatrixFSym &source1,
const TMatrixFSym &source2);
2074 template TMatrixFSym
operator+ <
Float_t>(
const TMatrixFSym &source1,
Float_t val);
2075 template TMatrixFSym
operator+ <
Float_t>(
Float_t val ,
const TMatrixFSym &source2);
2076 template TMatrixFSym
operator- <
Float_t>(
const TMatrixFSym &source1,
const TMatrixFSym &source2);
2077 template TMatrixFSym
operator- <
Float_t>(
const TMatrixFSym &source1,
Float_t val);
2078 template TMatrixFSym
operator- <
Float_t>(
Float_t val ,
const TMatrixFSym &source2);
2079 template TMatrixFSym
operator* <
Float_t>(
const TMatrixFSym &source,
Float_t val );
2080 template TMatrixFSym
operator* <
Float_t>(
Float_t val,
const TMatrixFSym &source );
2081 template TMatrixFSym
operator&& <
Float_t>(
const TMatrixFSym &source1,
const TMatrixFSym &source2);
2082 template TMatrixFSym
operator|| <
Float_t>(
const TMatrixFSym &source1,
const TMatrixFSym &source2);
2083 template TMatrixFSym
operator> <
Float_t>(
const TMatrixFSym &source1,
const TMatrixFSym &source2);
2084 template TMatrixFSym
operator>= <
Float_t>(
const TMatrixFSym &source1,
const TMatrixFSym &source2);
2085 template TMatrixFSym
operator<= <Float_t>(
const TMatrixFSym &source1,
const TMatrixFSym &source2);
2086 template TMatrixFSym
operator< <Float_t>(
const TMatrixFSym &source1,
const TMatrixFSym &source2);
2088 template TMatrixFSym &
Add <Float_t>(TMatrixFSym &target,
Float_t scalar,
const TMatrixFSym &source);
2092 #ifndef ROOT_TMatrixDSymfwd
2096 template class TMatrixTSym<Double_t>;
2098 template Bool_t operator== <
Double_t>(
const TMatrixDSym &source1,
const TMatrixDSym &source2);
2099 template TMatrixDSym
operator+ <
Double_t>(
const TMatrixDSym &source1,
const TMatrixDSym &source2);
2100 template TMatrixDSym
operator+ <
Double_t>(
const TMatrixDSym &source1,
Double_t val);
2101 template TMatrixDSym
operator+ <
Double_t>(
Double_t val ,
const TMatrixDSym &source2);
2102 template TMatrixDSym
operator- <
Double_t>(
const TMatrixDSym &source1,
const TMatrixDSym &source2);
2103 template TMatrixDSym
operator- <
Double_t>(
const TMatrixDSym &source1,
Double_t val);
2104 template TMatrixDSym
operator- <
Double_t>(
Double_t val ,
const TMatrixDSym &source2);
2105 template TMatrixDSym
operator* <
Double_t>(
const TMatrixDSym &source,
Double_t val );
2106 template TMatrixDSym
operator* <
Double_t>(
Double_t val,
const TMatrixDSym &source );
2107 template TMatrixDSym
operator&& <
Double_t>(
const TMatrixDSym &source1,
const TMatrixDSym &source2);
2108 template TMatrixDSym
operator|| <
Double_t>(
const TMatrixDSym &source1,
const TMatrixDSym &source2);
2109 template TMatrixDSym
operator> <
Double_t>(
const TMatrixDSym &source1,
const TMatrixDSym &source2);
2110 template TMatrixDSym
operator>= <
Double_t>(
const TMatrixDSym &source1,
const TMatrixDSym &source2);
2112 template TMatrixDSym
operator< <Double_t>(
const TMatrixDSym &source1,
const TMatrixDSym &source2);
TMatrixTSym< Element > & SimilarityT(const TMatrixT< Element > &n)
Calculate B^T * (*this) * B , final matrix will be (ncolsb x ncolsb) It is more efficient than applyi...
const TVectorD & GetEigenValues() const
virtual const Element * GetMatrixArray() const =0
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 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...
template TMatrixFSym operator<=< Float_t >(const TMatrixFSym &source1, const TMatrixFSym &source2)
int Invert(LASymMatrix &)
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
static Vc_ALWAYS_INLINE int_v min(const int_v &x, const int_v &y)
template TMatrixDSym operator<< Double_t >(const TMatrixDSym &source1, const TMatrixDSym &source2)
TMatrixTSym< Element > & ElementMult(TMatrixTSym< Element > &target, const TMatrixTSym< Element > &source)
Multiply target by the source, element-by-element.
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
TMatrixTSym< Element > & Use(Int_t row_lwb, Int_t row_upb, Element *data)
Small helper to encapsulate whether to return the value pointed to by the iterator or its address...
virtual TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *option="")
Copy array data to matrix .
TMatrixTSym< Element > & Add(TMatrixTSym< Element > &target, Element scalar, const TMatrixTSym< Element > &source)
Modify addition: target += scalar * source.
void Allocate(Int_t nrows, Int_t ncols, Int_t row_lwb=0, Int_t col_lwb=0, Int_t init=0, Int_t=-1)
Allocate new matrix.
TMatrixTSym< Element > & operator=(const TMatrixTSym< Element > &source)
templateClassImp(TMatrixTSym) template< class Element > TMatrixTSym< Element >
const TMatrixT< Element > EigenVectors(TVectorT< Element > &eigenValues) const
Return a matrix containing the eigen-vectors ordered by descending eigen-values.
TMatrixTSym< Element > & ElementDiv(TMatrixTSym< Element > &target, const TMatrixTSym< Element > &source)
Multiply target by the source, element-by-element.
TMatrixTSym< Element > operator>=(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
source1 >= source2
void ToUpper()
Change string to upper case.
Buffer base class used for serializing objects.
void TMult(const TMatrixT< Element > &a)
Create a matrix C such that C = A' * A.
Int_t GetNoElements() const
Expr< TransposeOp< SMatrix< T, D1, D2, R >, T, D1, D2 >, T, D2, D1, typename TranspPolicy< T, D1, D2, R >::RepType > Transpose(const SMatrix< T, D1, D2, R > &rhs)
Matrix Transpose B(i,j) = A(j,i) returning a matrix expression.
static Bool_t InvertLU(TMatrixD &a, Double_t tol, Double_t *det=0)
Calculate matrix inversion through in place forward/backward substitution.
Short_t Min(Short_t a, Short_t b)
virtual void Det(Double_t &d1, Double_t &d2)
Calculate determinant det = d1*TMath::Power(2.,d2)
TMatrixTSym< Element > & operator*=(Element val)
Multiply every element of the matrix with val.
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t nr_nonzeros=-1)=0
template TMatrixFSym operator<< Float_t >(const TMatrixFSym &source1, const TMatrixFSym &source2)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
TMatrixTSym< Element > & Rank1Update(const TVectorT< Element > &v, Element alpha=1.0)
Perform a rank 1 operation on the matrix: A += alpha * v * v^T.
double beta(double x, double y)
Calculates the beta function.
virtual Bool_t IsSymmetric() const
Check whether matrix is symmetric.
template TMatrixDSym & ElementDiv< Double_t >(TMatrixDSym &target, const TMatrixDSym &source)
TMatrixTSym< Element > & operator+=(Element val)
Add val to every element of the matrix.
virtual void Operation(Element &element) const =0
virtual Double_t Determinant() const
void AMultB(int n, int m, int k, const double *A, const double *B, double *C)
template TMatrixDSym & ElementMult< Double_t >(TMatrixDSym &target, const TMatrixDSym &source)
TMatrixTSym< Element > operator&&(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
Logical AND.
static double p2(double t, double a, double b, double c)
TObject & operator=(const TObject &rhs)
TObject assignment operator.
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 .
void TMult(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
Create a matrix C such that C = A' * B.
virtual TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *option="")
Copy array data to matrix .
virtual void FillIn(TMatrixTSym< Element > &m) const =0
void Error(const char *location, const char *msgfmt,...)
void Minus(const TMatrixTSym< Element > &a, const TMatrixTSym< Element > &b)
Symmetric matrix summation. Create a matrix C such that C = A + B.
virtual const Element * GetMatrixArray() const
TMatrixTSym< Element > operator>(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
source1 > source2
Int_t Memcpy_m(Element *newp, const Element *oldp, Int_t copySize, Int_t newSize, Int_t oldSize)
copy copySize doubles from *oldp to *newp .
Element * GetMatrixArray()
TMatrixTBase< Element > & Apply(const TElementActionT< Element > &action)
Apply action to each matrix element.
template TMatrixFSym & Add< Float_t >(TMatrixFSym &target, Float_t scalar, const TMatrixFSym &source)
template TMatrixDSym operator<=< Double_t >(const TMatrixDSym &source1, const TMatrixDSym &source2)
template TMatrixFSym & ElementDiv< Float_t >(TMatrixFSym &target, const TMatrixFSym &source)
virtual const Element * GetMatrixArray() const
static double p1(double t, double a, double b)
Bool_t operator==(const TMatrixTSym< Element > &m1, const TMatrixTSym< Element > &m2)
Check to see if two matrices are identical.
R__EXTERN Int_t gMatrixCheck
virtual TMatrixTBase< Element > & Randomize(Element alpha, Element beta, Double_t &seed)
randomize matrix element values but keep matrix symmetric
virtual void ReadFastArray(Bool_t *b, Int_t n)=0
virtual void WriteFastArray(const Bool_t *b, Int_t n)=0
TMatrixTSym< Element > & Transpose(const TMatrixTSym< Element > &source)
Transpose a matrix.
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
void Plus(const TMatrixTSym< Element > &a, const TMatrixTSym< Element > &b)
Symmetric matrix summation. Create a matrix C such that C = A + B.
TMatrixTSym< Element > & operator-=(Element val)
Subtract val from every element of the matrix.
const TMatrixD & GetEigenVectors() const
virtual TMatrixTSym< Element > & RandomizePD(Element alpha, Element beta, Double_t &seed)
randomize matrix element values but keep matrix symmetric positive definite
TMatrixTSym< Element > & InvertFast(Double_t *det=0)
Invert the matrix and calculate its determinant.
TMatrixTSym< Element > & SetSub(Int_t row_lwb, const TMatrixTBase< Element > &source)
Insert matrix source starting at [row_lwb][row_lwb], thereby overwriting the part [row_lwb...
TMatrixTSym< Element > & Similarity(const TMatrixT< Element > &n)
Calculate B * (*this) * B^T , final matrix will be (nrowsb x nrowsb) This is a similarity transform w...
virtual const Int_t * GetRowIndexArray() const
void Delete_m(Int_t size, Element *&)
delete data pointer m, if it was assigned on the heap
virtual void Operation(Element &element) const =0
TCppObject_t Allocate(TCppType_t type)
TMatrixT< Element > & Use(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, Element *data)
Use the array data to fill the matrix ([row_lwb..row_upb] x [col_lwb..col_upb])
TMatrixTSym< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
TMatrixTSym< Element > operator||(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
Logical Or.
Int_t GetNoElements() const
virtual const Int_t * GetColIndexArray() const
TMatrixTSym< Element > operator+(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
Short_t Max(Short_t a, Short_t b)
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
template TMatrixDSym & Add< Double_t >(TMatrixDSym &target, Double_t scalar, const TMatrixDSym &source)
Element * New_m(Int_t size)
return data pointer .
TMatrixTSym< Element > operator*(const TMatrixTSym< Element > &source1, Element val)
template TMatrixFSym & ElementMult< Float_t >(TMatrixFSym &target, const TMatrixFSym &source)
Element * fElements
data container
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
TMatrixTSym< Element > operator-(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
TMatrixTSym< Element > & GetSub(Int_t row_lwb, Int_t row_upb, TMatrixTSym< Element > &target, Option_t *option="S") const
Get submatrix [row_lwb..row_upb][row_lwb..row_upb]; The indexing range of the returned matrix depends...
virtual const Int_t * GetColIndexArray() const =0
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
virtual const Int_t * GetRowIndexArray() const =0