41template<
class Element>
49template<
class Element>
52 const Int_t no_rows = row_upb-row_lwb+1;
53 Allocate(no_rows,no_rows,row_lwb,row_lwb,1);
65template<
class Element>
69 SetMatrixArray(elements,option);
70 if (!this->IsSymmetric()) {
71 Error(
"TMatrixTSym(Int_t,Element*,Option_t*)",
"matrix not symmetric");
78template<
class Element>
81 const Int_t no_rows = row_upb-row_lwb+1;
82 Allocate(no_rows,no_rows,row_lwb,row_lwb);
83 SetMatrixArray(elements,option);
84 if (!this->IsSymmetric()) {
85 Error(
"TMatrixTSym(Int_t,Int_t,Element*,Option_t*)",
"matrix not symmetric");
91template<
class Element>
104template<
class Element>
134 const Element oldTol = this->SetTol(std::numeric_limits<Element>::min());
136 this->SetTol(oldTol);
146 Error(
"TMatrixTSym(EMatrixCreatorOp1,const TMatrixTSym)",
147 "operation %d not yet implemented", op);
153template<
class Element>
165 Error(
"TMatrixTSym(EMatrixCreatorOp1,const TMatrixT)",
166 "operation %d not yet implemented", op);
172template<
class Element>
181 Allocate(
a.GetNcols(),
a.GetNcols(),
a.GetColLwb(),
a.GetColLwb(),1);
188 Allocate(
a.GetNcols(),
a.GetNcols(),
a.GetColLwb(),
a.GetColLwb(),1);
194 Error(
"TMatrixTSym(EMatrixCreatorOp2)",
"operation %d not yet implemented", op);
200template<
class Element>
206 lazy_constructor.
FillIn(*
this);
207 if (!this->IsSymmetric()) {
208 Error(
"TMatrixTSym(TMatrixTSymLazy)",
"matrix not symmetric");
215template<
class Element>
219 if (size > this->kSizeMax)
229template<
class Element>
232 if (size == 0)
return 0;
234 if ( size <= this->kSizeMax )
237 Element *heap =
new Element[size];
247template<
class Element>
251 if (copySize == 0 || oldp == newp)
254 if ( newSize <= this->kSizeMax && oldSize <= this->kSizeMax ) {
257 for (
Int_t i = copySize-1; i >= 0; i--)
260 for (
Int_t i = 0; i < copySize; i++)
265 memcpy(newp,oldp,copySize*
sizeof(Element));
274template<
class Element>
278 this->fIsOwner =
kTRUE;
287 if (no_rows < 0 || no_cols < 0)
289 Error(
"Allocate",
"no_rows=%d no_cols=%d",no_rows,no_cols);
295 this->fNrows = no_rows;
296 this->fNcols = no_cols;
297 this->fRowLwb = row_lwb;
298 this->fColLwb = col_lwb;
299 this->fNelems = this->fNrows*this->fNcols;
301 if (this->fNelems > 0) {
302 fElements = New_m(this->fNelems);
304 memset(fElements,0,this->fNelems*
sizeof(Element));
312template<
class Element>
317 Error(
"Plus",
"matrices not compatible");
321 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
322 Error(
"Plus",
"this->GetMatrixArray() == a.GetMatrixArray()");
326 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
327 Error(
"Plus",
"this->GetMatrixArray() == b.GetMatrixArray()");
332 const Element * ap =
a.GetMatrixArray();
333 const Element * bp =
b.GetMatrixArray();
334 Element * cp = this->GetMatrixArray();
335 const Element *
const cp_last = cp+this->fNelems;
337 while (cp < cp_last) {
346template<
class Element>
351 Error(
"Minus",
"matrices not compatible");
355 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
356 Error(
"Minus",
"this->GetMatrixArray() == a.GetMatrixArray()");
360 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
361 Error(
"Minus",
"this->GetMatrixArray() == b.GetMatrixArray()");
366 const Element * ap =
a.GetMatrixArray();
367 const Element * bp =
b.GetMatrixArray();
368 Element * cp = this->GetMatrixArray();
369 const Element *
const cp_last = cp+this->fNelems;
371 while (cp < cp_last) {
381template<
class Element>
387 const Element *ap =
a.GetMatrixArray();
388 Element *cp = this->GetMatrixArray();
389 if (
typeid(Element) ==
typeid(
Double_t))
390 cblas_dgemm (CblasRowMajor,CblasTrans,CblasNoTrans,this->fNrows,this->fNcols,
a.GetNrows(),
391 1.0,ap,
a.GetNcols(),ap,
a.GetNcols(),1.0,cp,this->fNcols);
392 else if (
typeid(Element) !=
typeid(
Float_t))
393 cblas_sgemm (CblasRowMajor,CblasTrans,CblasNoTrans,fNrows,fNcols,
a.GetNrows(),
394 1.0,ap,
a.GetNcols(),ap,
a.GetNcols(),1.0,cp,fNcols);
396 Error(
"TMult",
"type %s not implemented in BLAS library",
typeid(Element));
398 const Int_t nb =
a.GetNoElements();
399 const Int_t ncolsa =
a.GetNcols();
400 const Int_t ncolsb = ncolsa;
401 const Element *
const ap =
a.GetMatrixArray();
402 const Element *
const bp = ap;
403 Element * cp = this->GetMatrixArray();
405 const Element *acp0 = ap;
406 while (acp0 < ap+
a.GetNcols()) {
407 for (
const Element *bcp = bp; bcp < bp+ncolsb; ) {
408 const Element *acp = acp0;
410 while (bcp < bp+nb) {
421 R__ASSERT(cp == this->GetMatrixArray()+this->fNelems && acp0 == ap+ncolsa);
429template<
class Element>
435 const Element *ap =
a.GetMatrixArray();
436 Element *cp = this->GetMatrixArray();
437 if (
typeid(Element) ==
typeid(
Double_t))
438 cblas_dsymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
439 ap,
a.GetNcols(),ap,
a.GetNcols(),0.0,cp,this->fNcols);
440 else if (
typeid(Element) !=
typeid(
Float_t))
441 cblas_ssymm (CblasRowMajor,CblasLeft,CblasUpper,fNrows,fNcols,1.0,
442 ap1,
a.GetNcols(),ap,
a.GetNcols(),0.0,cp,fNcols);
444 Error(
"TMult",
"type %s not implemented in BLAS library",
typeid(Element));
446 const Int_t nb =
a.GetNoElements();
447 const Int_t ncolsa =
a.GetNcols();
448 const Int_t ncolsb = ncolsa;
449 const Element *
const ap =
a.GetMatrixArray();
450 const Element *
const bp = ap;
451 Element * cp = this->GetMatrixArray();
453 const Element *acp0 = ap;
454 while (acp0 < ap+
a.GetNcols()) {
455 for (
const Element *bcp = bp; bcp < bp+ncolsb; ) {
456 const Element *acp = acp0;
458 while (bcp < bp+nb) {
469 R__ASSERT(cp == this->GetMatrixArray()+this->fNelems && acp0 == ap+ncolsa);
475template<
class Element>
480 Error(
"Use",
"row_upb=%d < row_lwb=%d",row_upb,row_lwb);
485 this->fNrows = row_upb-row_lwb+1;
486 this->fNcols = this->fNrows;
487 this->fRowLwb = row_lwb;
488 this->fColLwb = row_lwb;
489 this->fNelems = this->fNrows*this->fNcols;
503template<
class Element>
509 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
510 Error(
"GetSub",
"row_lwb out of bounds");
513 if (row_upb < this->fRowLwb || row_upb > this->fRowLwb+this->fNrows-1) {
514 Error(
"GetSub",
"row_upb out of bounds");
517 if (row_upb < row_lwb) {
518 Error(
"GetSub",
"row_upb < row_lwb");
531 row_upb_sub = row_upb-row_lwb;
533 row_lwb_sub = row_lwb;
534 row_upb_sub = row_upb;
537 target.
ResizeTo(row_lwb_sub,row_upb_sub,row_lwb_sub,row_upb_sub);
538 const Int_t nrows_sub = row_upb_sub-row_lwb_sub+1;
541 for (
Int_t irow = 0; irow < nrows_sub; irow++) {
542 for (
Int_t icol = 0; icol < nrows_sub; icol++) {
543 target(irow+row_lwb_sub,icol+row_lwb_sub) = (*this)(row_lwb+irow,row_lwb+icol);
547 const Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNrows+(row_lwb-this->fRowLwb);
550 for (
Int_t irow = 0; irow < nrows_sub; irow++) {
551 const Element *ap_sub = ap;
552 for (
Int_t icol = 0; icol < nrows_sub; icol++) {
569template<
class Element>
575 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
576 Error(
"GetSub",
"row_lwb out of bounds");
579 if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
580 Error(
"GetSub",
"col_lwb out of bounds");
583 if (row_upb < this->fRowLwb || row_upb > this->fRowLwb+this->fNrows-1) {
584 Error(
"GetSub",
"row_upb out of bounds");
587 if (col_upb < this->fColLwb || col_upb > this->fColLwb+this->fNcols-1) {
588 Error(
"GetSub",
"col_upb out of bounds");
591 if (row_upb < row_lwb || col_upb < col_lwb) {
592 Error(
"GetSub",
"row_upb < row_lwb || col_upb < col_lwb");
601 const Int_t row_lwb_sub = (shift) ? 0 : row_lwb;
602 const Int_t row_upb_sub = (shift) ? row_upb-row_lwb : row_upb;
603 const Int_t col_lwb_sub = (shift) ? 0 : col_lwb;
604 const Int_t col_upb_sub = (shift) ? col_upb-col_lwb : col_upb;
606 target.
ResizeTo(row_lwb_sub,row_upb_sub,col_lwb_sub,col_upb_sub);
607 const Int_t nrows_sub = row_upb_sub-row_lwb_sub+1;
608 const Int_t ncols_sub = col_upb_sub-col_lwb_sub+1;
611 for (
Int_t irow = 0; irow < nrows_sub; irow++) {
612 for (
Int_t icol = 0; icol < ncols_sub; icol++) {
613 target(irow+row_lwb_sub,icol+col_lwb_sub) = (*this)(row_lwb+irow,col_lwb+icol);
617 const Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNcols+(col_lwb-this->fColLwb);
620 for (
Int_t irow = 0; irow < nrows_sub; irow++) {
621 const Element *ap_sub = ap;
622 for (
Int_t icol = 0; icol < ncols_sub; icol++) {
636template<
class Element>
644 Error(
"SetSub",
"source matrix is not symmetric");
647 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
648 Error(
"SetSub",
"row_lwb outof bounds");
651 if (row_lwb+source.
GetNrows() > this->fRowLwb+this->fNrows) {
652 Error(
"SetSub",
"source matrix too large");
661 for (
Int_t irow = 0; irow < nRows_source; irow++) {
662 for (
Int_t icol = 0; icol < nRows_source; icol++) {
663 (*this)(row_lwb+irow,row_lwb+icol) = source(rowlwb_s+irow,rowlwb_s+icol);
668 Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNrows+(row_lwb-this->fRowLwb);
670 for (
Int_t irow = 0; irow < nRows_source; irow++) {
671 Element *ap_sub = ap;
672 for (
Int_t icol = 0; icol < nRows_source; icol++) {
686template<
class Element>
693 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
694 Error(
"SetSub",
"row_lwb out of bounds");
697 if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
698 Error(
"SetSub",
"col_lwb out of bounds");
702 if (row_lwb+source.
GetNrows() > this->fRowLwb+this->fNrows || col_lwb+source.
GetNcols() > this->fRowLwb+this->fNrows) {
703 Error(
"SetSub",
"source matrix too large");
706 if (col_lwb+source.
GetNcols() > this->fRowLwb+this->fNrows || row_lwb+source.
GetNrows() > this->fRowLwb+this->fNrows) {
707 Error(
"SetSub",
"source matrix too large");
717 if (row_lwb >= col_lwb) {
720 for (irow = 0; irow < nRows_source; irow++) {
721 for (
Int_t icol = 0; col_lwb+icol <= row_lwb+irow &&
722 icol < nCols_source; icol++) {
723 (*this)(row_lwb+irow,col_lwb+icol) = source(irow+rowlwb_s,icol+collwb_s);
728 for (irow = 0; irow < nCols_source; irow++) {
729 for (
Int_t icol = nRows_source-1; row_lwb+icol > irow+col_lwb &&
731 (*this)(col_lwb+irow,row_lwb+icol) = source(icol+rowlwb_s,irow+collwb_s);
743template<
class Element>
747 if (!this->IsSymmetric()) {
748 Error(
"SetMatrixArray",
"Matrix is not symmetric after Set");
756template<
class Element>
759 if (row_shift != col_shift) {
760 Error(
"Shift",
"row_shift != col_shift");
771template<
class Element>
775 if (!this->fIsOwner) {
776 Error(
"ResizeTo(Int_t,Int_t)",
"Not owner of data array,cannot resize");
780 if (nrows != ncols) {
781 Error(
"ResizeTo(Int_t,Int_t)",
"nrows != ncols");
785 if (this->fNelems > 0) {
786 if (this->fNrows == nrows && this->fNcols == ncols)
788 else if (nrows == 0 || ncols == 0) {
789 this->fNrows = nrows; this->fNcols = ncols;
794 Element *elements_old = GetMatrixArray();
795 const Int_t nelems_old = this->fNelems;
796 const Int_t nrows_old = this->fNrows;
797 const Int_t ncols_old = this->fNcols;
802 Element *elements_new = GetMatrixArray();
805 if (this->fNelems > this->kSizeMax || nelems_old > this->kSizeMax)
806 memset(elements_new,0,this->fNelems*
sizeof(Element));
807 else if (this->fNelems > nelems_old)
808 memset(elements_new+nelems_old,0,(this->fNelems-nelems_old)*
sizeof(Element));
814 const Int_t nelems_new = this->fNelems;
815 if (ncols_old < this->fNcols) {
816 for (
Int_t i = nrows_copy-1; i >= 0; i--) {
817 Memcpy_m(elements_new+i*this->fNcols,elements_old+i*ncols_old,ncols_copy,
818 nelems_new,nelems_old);
819 if (this->fNelems <= this->kSizeMax && nelems_old <= this->kSizeMax)
820 memset(elements_new+i*this->fNcols+ncols_copy,0,(this->fNcols-ncols_copy)*
sizeof(Element));
823 for (
Int_t i = 0; i < nrows_copy; i++)
824 Memcpy_m(elements_new+i*this->fNcols,elements_old+i*ncols_old,ncols_copy,
825 nelems_new,nelems_old);
828 Delete_m(nelems_old,elements_old);
841template<
class Element>
846 if (!this->fIsOwner) {
847 Error(
"ResizeTo(Int_t,Int_t,Int_t,Int_t)",
"Not owner of data array,cannot resize");
851 if (row_lwb != col_lwb) {
852 Error(
"ResizeTo(Int_t,Int_t,Int_t,Int_t)",
"row_lwb != col_lwb");
855 if (row_upb != col_upb) {
856 Error(
"ResizeTo(Int_t,Int_t,Int_t,Int_t)",
"row_upb != col_upb");
860 const Int_t new_nrows = row_upb-row_lwb+1;
861 const Int_t new_ncols = col_upb-col_lwb+1;
863 if (this->fNelems > 0) {
865 if (this->fNrows == new_nrows && this->fNcols == new_ncols &&
866 this->fRowLwb == row_lwb && this->fColLwb == col_lwb)
868 else if (new_nrows == 0 || new_ncols == 0) {
869 this->fNrows = new_nrows; this->fNcols = new_ncols;
870 this->fRowLwb = row_lwb; this->fColLwb = col_lwb;
875 Element *elements_old = GetMatrixArray();
876 const Int_t nelems_old = this->fNelems;
877 const Int_t nrows_old = this->fNrows;
878 const Int_t ncols_old = this->fNcols;
879 const Int_t rowLwb_old = this->fRowLwb;
880 const Int_t colLwb_old = this->fColLwb;
882 Allocate(new_nrows,new_ncols,row_lwb,col_lwb);
885 Element *elements_new = GetMatrixArray();
888 if (this->fNelems > this->kSizeMax || nelems_old > this->kSizeMax)
889 memset(elements_new,0,this->fNelems*
sizeof(Element));
890 else if (this->fNelems > nelems_old)
891 memset(elements_new+nelems_old,0,(this->fNelems-nelems_old)*
sizeof(Element));
896 const Int_t rowUpb_copy =
TMath::Min(this->fRowLwb+this->fNrows-1,rowLwb_old+nrows_old-1);
897 const Int_t colUpb_copy =
TMath::Min(this->fColLwb+this->fNcols-1,colLwb_old+ncols_old-1);
899 const Int_t nrows_copy = rowUpb_copy-rowLwb_copy+1;
900 const Int_t ncols_copy = colUpb_copy-colLwb_copy+1;
902 if (nrows_copy > 0 && ncols_copy > 0) {
903 const Int_t colOldOff = colLwb_copy-colLwb_old;
904 const Int_t colNewOff = colLwb_copy-this->fColLwb;
905 if (ncols_old < this->fNcols) {
906 for (
Int_t i = nrows_copy-1; i >= 0; i--) {
907 const Int_t iRowOld = rowLwb_copy+i-rowLwb_old;
908 const Int_t iRowNew = rowLwb_copy+i-this->fRowLwb;
909 Memcpy_m(elements_new+iRowNew*this->fNcols+colNewOff,
910 elements_old+iRowOld*ncols_old+colOldOff,ncols_copy,this->fNelems,nelems_old);
911 if (this->fNelems <= this->kSizeMax && nelems_old <= this->kSizeMax)
912 memset(elements_new+iRowNew*this->fNcols+colNewOff+ncols_copy,0,
913 (this->fNcols-ncols_copy)*
sizeof(Element));
916 for (
Int_t i = 0; i < nrows_copy; i++) {
917 const Int_t iRowOld = rowLwb_copy+i-rowLwb_old;
918 const Int_t iRowNew = rowLwb_copy+i-this->fRowLwb;
919 Memcpy_m(elements_new+iRowNew*this->fNcols+colNewOff,
920 elements_old+iRowOld*ncols_old+colOldOff,ncols_copy,this->fNelems,nelems_old);
925 Delete_m(nelems_old,elements_old);
927 Allocate(new_nrows,new_ncols,row_lwb,col_lwb,1);
935template<
class Element>
947template<
class Element>
961template<
class Element>
968 Element *p2 = this->GetMatrixArray();
969 for (
Int_t i = 0; i < this->GetNoElements(); i++)
979template<
class Element>
988 Element *pM = this->GetMatrixArray();
990 Error(
"InvertFast",
"matrix is singular");
1000 TMatrixTSymCramerInv::Inv2x2<Element>(*
this,det);
1005 TMatrixTSymCramerInv::Inv3x3<Element>(*
this,det);
1010 TMatrixTSymCramerInv::Inv4x4<Element>(*
this,det);
1015 TMatrixTSymCramerInv::Inv5x5<Element>(*
this,det);
1020 TMatrixTSymCramerInv::Inv6x6<Element>(*
this,det);
1029 Element *p2 = this->GetMatrixArray();
1030 for (
Int_t i = 0; i < this->GetNoElements(); i++)
1041template<
class Element>
1050 Error(
"Transpose",
"matrix has wrong shape");
1063template<
class Element>
1069 if (
v.GetNoElements() < this->fNrows) {
1070 Error(
"Rank1Update",
"vector too short");
1075 const Element *
const pv =
v.GetMatrixArray();
1076 Element *trp = this->GetMatrixArray();
1078 for (
Int_t i = 0; i < this->fNrows; i++) {
1080 tcp += i*this->fNcols;
1081 const Element tmp = alpha*pv[i];
1082 for (
Int_t j = i; j < this->fNcols; j++) {
1083 if (j > i) *tcp += tmp*pv[j];
1084 *trp++ += tmp*pv[j];
1085 tcp += this->fNcols;
1087 tcp -= this->fNelems-1;
1099template<
class Element>
1105 if (this->fNcols !=
b.GetNcols() || this->fColLwb !=
b.GetColLwb()) {
1106 Error(
"Similarity(const TMatrixT &)",
"matrices incompatible");
1111 const Int_t ncolsa = this->fNcols;
1112 const Int_t nb =
b.GetNoElements();
1113 const Int_t nrowsb =
b.GetNrows();
1114 const Int_t ncolsb =
b.GetNcols();
1116 const Element *
const bp =
b.GetMatrixArray();
1118 Element work[kWorkMax];
1120 Element *bap = work;
1121 if (nrowsb*ncolsa > kWorkMax) {
1122 isAllocated =
kTRUE;
1123 bap =
new Element[nrowsb*ncolsa];
1126 AMultB(bp,nb,ncolsb,this->fElements,this->fNelems,this->fNcols,bap);
1128 if (nrowsb != this->fNrows)
1129 this->ResizeTo(nrowsb,nrowsb);
1132 Element *cp = this->GetMatrixArray();
1133 if (
typeid(Element) ==
typeid(
Double_t))
1134 cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasTrans,this->fNrows,this->fNcols,ba.GetNcols(),
1135 1.0,bap,ba.GetNcols(),bp,
b.GetNcols(),1.0,cp,this->fNcols);
1136 else if (
typeid(Element) !=
typeid(
Float_t))
1137 cblas_sgemm (CblasRowMajor,CblasNoTrans,CblasTrans,this->fNrows,this->fNcols,ba.GetNcols(),
1138 1.0,bap,ba.GetNcols(),bp,
b.GetNcols(),1.0,cp,this->fNcols);
1140 Error(
"Similarity",
"type %s not implemented in BLAS library",
typeid(Element));
1142 const Int_t nba = nrowsb*ncolsa;
1143 const Int_t ncolsba = ncolsa;
1144 const Element * bi1p = bp;
1145 Element * cp = this->GetMatrixArray();
1146 Element *
const cp0 = cp;
1149 const Element *barp0 = bap;
1150 while (barp0 < bap+nba) {
1151 const Element *brp0 = bi1p;
1152 while (brp0 < bp+nb) {
1153 const Element *barp = barp0;
1154 const Element *brp = brp0;
1156 while (brp < brp0+ncolsb)
1157 cij += *barp++ * *brp++;
1166 R__ASSERT(cp == cp0+this->fNelems+ishift && barp0 == bap+nba);
1169 for (
Int_t irow = 0; irow < this->fNrows; irow++) {
1170 const Int_t rowOff1 = irow*this->fNrows;
1171 for (
Int_t icol = 0; icol < irow; icol++) {
1172 const Int_t rowOff2 = icol*this->fNrows;
1173 cp[rowOff1+icol] = cp[rowOff2+irow];
1190template<
class Element>
1196 if (this->fNcols !=
b.GetNcols() || this->fColLwb !=
b.GetColLwb()) {
1197 Error(
"Similarity(const TMatrixTSym &)",
"matrices incompatible");
1203 const Int_t nrowsb =
b.GetNrows();
1204 const Int_t ncolsa = this->GetNcols();
1206 Element work[kWorkMax];
1208 Element *abtp = work;
1209 if (this->fNcols > kWorkMax) {
1210 isAllocated =
kTRUE;
1211 abtp =
new Element[this->fNcols];
1216 const Element *bp =
b.GetMatrixArray();
1217 Element *cp = this->GetMatrixArray();
1218 if (
typeid(Element) ==
typeid(
Double_t))
1219 cblas_dsymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
1220 bp,
b.GetNcols(),abtp,abt.
GetNcols(),0.0,cp,this->fNcols);
1221 else if (
typeid(Element) !=
typeid(
Float_t))
1222 cblas_ssymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
1223 bp,
b.GetNcols(),abtp,abt.
GetNcols(),0.0,cp,this->fNcols);
1225 Error(
"Similarity",
"type %s not implemented in BLAS library",
typeid(Element));
1230 const Int_t ncolsa = this->GetNcols();
1231 const Int_t nb =
b.GetNoElements();
1232 const Int_t nrowsb =
b.GetNrows();
1233 const Int_t ncolsb =
b.GetNcols();
1235 const Element *
const bp =
b.GetMatrixArray();
1237 Element work[kWorkMax];
1239 Element *bap = work;
1240 if (nrowsb*ncolsa > kWorkMax) {
1241 isAllocated =
kTRUE;
1242 bap =
new Element[nrowsb*ncolsa];
1245 AMultB(bp,nb,ncolsb,this->fElements,this->fNelems,this->fNcols,bap);
1247 const Int_t nba = nrowsb*ncolsa;
1248 const Int_t ncolsba = ncolsa;
1249 const Element * bi1p = bp;
1250 Element * cp = this->GetMatrixArray();
1251 Element *
const cp0 = cp;
1254 const Element *barp0 = bap;
1255 while (barp0 < bap+nba) {
1256 const Element *brp0 = bi1p;
1257 while (brp0 < bp+nb) {
1258 const Element *barp = barp0;
1259 const Element *brp = brp0;
1261 while (brp < brp0+ncolsb)
1262 cij += *barp++ * *brp++;
1271 R__ASSERT(cp == cp0+this->fNelems+ishift && barp0 == bap+nba);
1274 for (
Int_t irow = 0; irow < this->fNrows; irow++) {
1275 const Int_t rowOff1 = irow*this->fNrows;
1276 for (
Int_t icol = 0; icol < irow; icol++) {
1277 const Int_t rowOff2 = icol*this->fNrows;
1278 cp[rowOff1+icol] = cp[rowOff2+irow];
1292template<
class Element>
1298 if (this->fNcols !=
v.GetNrows() || this->fColLwb !=
v.GetLwb()) {
1299 Error(
"Similarity(const TVectorT &)",
"vector and matrix incompatible");
1304 const Element *mp = this->GetMatrixArray();
1305 const Element *vp =
v.GetMatrixArray();
1308 const Element *
const vp_first = vp;
1309 const Element *
const vp_last = vp+
v.GetNrows();
1310 while (vp < vp_last) {
1312 for (
const Element *sp = vp_first; sp < vp_last; )
1313 sum2 += *mp++ * *sp++;
1314 sum1 += sum2 * *vp++;
1317 R__ASSERT(mp == this->GetMatrixArray()+this->GetNoElements());
1327template<
class Element>
1333 if (this->fNrows !=
b.GetNrows() || this->fRowLwb !=
b.GetRowLwb()) {
1334 Error(
"SimilarityT(const TMatrixT &)",
"matrices incompatible");
1339 const Int_t ncolsb =
b.GetNcols();
1340 const Int_t ncolsa = this->GetNcols();
1342 Element work[kWorkMax];
1344 Element *btap = work;
1345 if (ncolsb*ncolsa > kWorkMax) {
1346 isAllocated =
kTRUE;
1347 btap =
new Element[ncolsb*ncolsa];
1353 if (ncolsb != this->fNcols)
1354 this->ResizeTo(ncolsb,ncolsb);
1357 const Element *bp =
b.GetMatrixArray();
1358 Element *cp = this->GetMatrixArray();
1359 if (
typeid(Element) ==
typeid(
Double_t))
1360 cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasNoTrans,this->fNrows,this->fNcols,bta.
GetNcols(),
1361 1.0,btap,bta.
GetNcols(),bp,
b.GetNcols(),1.0,cp,this->fNcols);
1362 else if (
typeid(Element) !=
typeid(
Float_t))
1363 cblas_sgemm (CblasRowMajor,CblasNoTrans,CblasNoTrans,this->fNrows,this->fNcols,bta.
GetNcols(),
1364 1.0,btap,bta.
GetNcols(),bp,
b.GetNcols(),1.0,cp,this->fNcols);
1366 Error(
"similarityT",
"type %s not implemented in BLAS library",
typeid(Element));
1369 const Int_t nb =
b.GetNoElements();
1371 const Element *
const bp =
b.GetMatrixArray();
1372 Element * cp = this->GetMatrixArray();
1373 Element *
const cp0 = cp;
1376 const Element *btarp0 = btap;
1377 const Element *bcp0 = bp;
1378 while (btarp0 < btap+nbta) {
1379 for (
const Element *bcp = bcp0; bcp < bp+ncolsb; ) {
1380 const Element *btarp = btarp0;
1382 while (bcp < bp+nb) {
1383 cij += *btarp++ * *bcp;
1394 R__ASSERT(cp == cp0+this->fNelems+ishift && btarp0 == btap+nbta);
1397 for (
Int_t irow = 0; irow < this->fNrows; irow++) {
1398 const Int_t rowOff1 = irow*this->fNrows;
1399 for (
Int_t icol = 0; icol < irow; icol++) {
1400 const Int_t rowOff2 = icol*this->fNrows;
1401 cp[rowOff1+icol] = cp[rowOff2+irow];
1414template<
class Element>
1418 Error(
"operator=",
"matrices not compatible");
1424 memcpy(this->GetMatrixArray(),source.
fElements,this->fNelems*
sizeof(Element));
1431template<
class Element>
1436 if (lazy_constructor.
fRowUpb != this->GetRowUpb() ||
1437 lazy_constructor.
fRowLwb != this->GetRowLwb()) {
1438 Error(
"operator=(const TMatrixTSymLazy&)",
"matrix is incompatible with "
1439 "the assigned Lazy matrix");
1443 lazy_constructor.
FillIn(*
this);
1450template<
class Element>
1455 Element *ep = fElements;
1456 const Element *
const ep_last = ep+this->fNelems;
1457 while (ep < ep_last)
1466template<
class Element>
1471 Element *ep = fElements;
1472 const Element *
const ep_last = ep+this->fNelems;
1473 while (ep < ep_last)
1482template<
class Element>
1487 Element *ep = fElements;
1488 const Element *
const ep_last = ep+this->fNelems;
1489 while (ep < ep_last)
1498template<
class Element>
1503 Element *ep = fElements;
1504 const Element *
const ep_last = ep+this->fNelems;
1505 while (ep < ep_last)
1514template<
class Element>
1518 Error(
"operator+=",
"matrices not compatible");
1523 Element *tp = this->GetMatrixArray();
1524 const Element *
const tp_last = tp+this->fNelems;
1525 while (tp < tp_last)
1534template<
class Element>
1538 Error(
"operator-=",
"matrices not compatible");
1543 Element *tp = this->GetMatrixArray();
1544 const Element *
const tp_last = tp+this->fNelems;
1545 while (tp < tp_last)
1553template<
class Element>
1559 Element *trp = this->GetMatrixArray();
1561 for (
Int_t i = 0; i < this->fNrows; i++) {
1563 tcp += i*this->fNcols;
1564 for (
Int_t j = i; j < this->fNcols; j++) {
1566 if (j > i) *tcp = val;
1568 tcp += this->fNcols;
1570 tcp -= this->fNelems-1;
1580template<
class Element>
1586 Element *trp = this->GetMatrixArray();
1588 for (
Int_t i = 0; i < this->fNrows; i++) {
1589 action.
fI = i+this->fRowLwb;
1591 tcp += i*this->fNcols;
1592 for (
Int_t j = i; j < this->fNcols; j++) {
1593 action.
fJ = j+this->fColLwb;
1595 if (j > i) *tcp = val;
1597 tcp += this->fNcols;
1599 tcp -= this->fNelems-1;
1608template<
class Element>
1613 if (this->fNrows != this->fNcols || this->fRowLwb != this->fColLwb) {
1614 Error(
"Randomize(Element,Element,Element &",
"matrix should be square");
1619 const Element scale =
beta-alpha;
1620 const Element shift = alpha/scale;
1622 Element *ep = GetMatrixArray();
1623 for (
Int_t i = 0; i < this->fNrows; i++) {
1624 const Int_t off = i*this->fNcols;
1625 for (
Int_t j = 0; j <= i; j++) {
1626 ep[off+j] = scale*(
Drand(seed)+shift);
1628 ep[j*this->fNcols+i] = ep[off+j];
1639template<
class Element>
1644 if (this->fNrows != this->fNcols || this->fRowLwb != this->fColLwb) {
1645 Error(
"RandomizeSym(Element,Element,Element &",
"matrix should be square");
1650 const Element scale =
beta-alpha;
1651 const Element shift = alpha/scale;
1653 Element *ep = GetMatrixArray();
1655 for (i = 0; i < this->fNrows; i++) {
1656 const Int_t off = i*this->fNcols;
1657 for (
Int_t j = 0; j <= i; j++)
1658 ep[off+j] = scale*(
Drand(seed)+shift);
1661 for (i = this->fNrows-1; i >= 0; i--) {
1662 const Int_t off1 = i*this->fNcols;
1663 for (
Int_t j = i; j >= 0; j--) {
1664 const Int_t off2 = j*this->fNcols;
1665 ep[off1+j] *= ep[off2+j];
1666 for (
Int_t k = j-1; k >= 0; k--) {
1667 ep[off1+j] += ep[off1+k]*ep[off2+k];
1670 ep[off2+i] = ep[off1+j];
1681template<
class Element>
1686 eigenValues.
ResizeTo(this->fNrows);
1694template<
class Element>
1704template<
class Element>
1714template<
class Element>
1724template<
class Element>
1732template<
class Element>
1742template<
class Element>
1752template<
class Element>
1755 return Element(-1.0)*
operator-(source1,val);
1760template<
class Element>
1770template<
class Element>
1779template<
class Element>
1785 Error(
"operator&&(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1795 while (tp < tp_last)
1796 *tp++ = (*sp1++ != 0.0 && *sp2++ != 0.0);
1804template<
class Element>
1810 Error(
"operator||(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1820 while (tp < tp_last)
1821 *tp++ = (*sp1++ != 0.0 || *sp2++ != 0.0);
1829template<
class Element>
1835 Error(
"operator>(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1845 while (tp < tp_last) {
1846 *tp++ = (*sp1) > (*sp2); sp1++; sp2++;
1855template<
class Element>
1861 Error(
"operator>=(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1871 while (tp < tp_last) {
1872 *tp++ = (*sp1) >= (*sp2); sp1++; sp2++;
1881template<
class Element>
1887 Error(
"operator<=(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1897 while (tp < tp_last) {
1898 *tp++ = (*sp1) <= (*sp2); sp1++; sp2++;
1907template<
class Element>
1913 Error(
"operator<(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1923 while (tp < tp_last) {
1924 *tp++ = (*sp1) < (*sp2); sp1++; sp2++;
1933template<
class Element>
1937 ::Error(
"Add",
"matrices not compatible");
1947 for (
Int_t i = 0; i < nrows; i++) {
1951 for (
Int_t j = i; j < ncols; j++) {
1952 const Element tmp = scalar * *sp++;
1953 if (j > i) *tcp += tmp;
1966template<
class Element>
1970 ::Error(
"ElementMult",
"matrices not compatible");
1980 for (
Int_t i = 0; i < nrows; i++) {
1984 for (
Int_t j = i; j < ncols; j++) {
1985 if (j > i) *tcp *= *sp;
1998template<
class Element>
2002 ::Error(
"ElementDiv",
"matrices not compatible");
2012 for (
Int_t i = 0; i < nrows; i++) {
2016 for (
Int_t j = i; j < ncols; j++) {
2018 if (j > i) *tcp /= *sp;
2023 Error(
"ElementDiv",
"source (%d,%d) is zero",irow,icol);
2037template<
class Element>
2045 fElements =
new Element[this->fNelems];
2047 for (i = 0; i < this->fNrows; i++) {
2048 R__b.
ReadFastArray(fElements+i*this->fNcols+i,this->fNcols-i);
2051 for (i = 0; i < this->fNrows; i++) {
2052 for (
Int_t j = 0; j < i; j++) {
2053 fElements[i*this->fNcols+j] = fElements[j*this->fNrows+i];
2056 if (this->fNelems <= this->kSizeMax) {
2057 memcpy(fDataStack,fElements,this->fNelems*
sizeof(Element));
2058 delete [] fElements;
2059 fElements = fDataStack;
2064 for (
Int_t i = 0; i < this->fNrows; i++) {
#define templateClassImp(name)
void Error(const char *location, const char *msgfmt,...)
R__EXTERN Int_t gMatrixCheck
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 .
template TMatrixDSym & ElementMult< Double_t >(TMatrixDSym &target, const TMatrixDSym &source)
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, Element val)
TMatrixTSym< Element > operator<(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
source1 < source2
template TMatrixFSym & ElementMult< Float_t >(TMatrixFSym &target, const TMatrixFSym &source)
TMatrixTSym< Element > operator+(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
template TMatrixFSym & Add< Float_t >(TMatrixFSym &target, Float_t scalar, const TMatrixFSym &source)
template TMatrixFSym & ElementDiv< Float_t >(TMatrixFSym &target, const TMatrixFSym &source)
template TMatrixDSym & Add< Double_t >(TMatrixDSym &target, Double_t scalar, const TMatrixDSym &source)
TMatrixTSym< Element > operator-(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
template TMatrixFSym operator<=< Float_t >(const TMatrixFSym &source1, const TMatrixFSym &source2)
TMatrixTSym< Element > & Add(TMatrixTSym< Element > &target, Element scalar, const TMatrixTSym< Element > &source)
Modify addition: target += scalar * source.
TMatrixTSym< Element > & ElementMult(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
TMatrixTSym< Element > operator>=(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
source1 >= source2
TMatrixTSym< Element > operator<=(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
source1 <= source2
template TMatrixFSym operator<<Float_t >(const TMatrixFSym &source1, const TMatrixFSym &source2)
template TMatrixDSym & ElementDiv< Double_t >(TMatrixDSym &target, const TMatrixDSym &source)
TMatrixTSym< Element > operator||(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
Logical Or.
TMatrixTSym< Element > operator&&(const TMatrixTSym< Element > &source1, const TMatrixTSym< Element > &source2)
Logical AND.
template TMatrixDSym operator<<Double_t >(const TMatrixDSym &source1, const TMatrixDSym &source2)
Bool_t operator==(const TMatrixTSym< Element > &m1, const TMatrixTSym< Element > &m2)
Check to see if two matrices are identical.
template TMatrixDSym operator<=< Double_t >(const TMatrixDSym &source1, const TMatrixDSym &source2)
Double_t Drand(Double_t &ix)
Random number generator [0....1] with seed ix.
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 void ReadFastArray(Bool_t *b, Int_t n)=0
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
virtual void WriteFastArray(const Bool_t *b, Int_t n)=0
virtual void Det(Double_t &d1, Double_t &d2)
Calculate determinant det = d1*TMath::Power(2.,d2)
static Bool_t InvertLU(TMatrixD &a, Double_t tol, Double_t *det=0)
Calculate matrix inversion through in place forward/backward substitution.
virtual void Operation(Element &element) const =0
virtual void Operation(Element &element) const =0
const TVectorD & GetEigenValues() const
const TMatrixD & GetEigenVectors() const
virtual const Element * GetMatrixArray() const =0
virtual const Int_t * GetRowIndexArray() const =0
virtual const Int_t * GetColIndexArray() const =0
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t nr_nonzeros=-1)=0
Int_t GetNoElements() const
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 Bool_t IsSymmetric() const
Check whether matrix is symmetric.
virtual TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *option="")
Copy array data to matrix .
virtual void FillIn(TMatrixTSym< Element > &m) const =0
TMatrixTSym< Element > & operator*=(Element val)
Multiply every element of the matrix with val.
virtual const Element * GetMatrixArray() const
TMatrixTSym< Element > & operator+=(Element val)
Add val to every element of the matrix.
virtual Double_t Determinant() const
virtual TMatrixTSym< Element > & RandomizePD(Element alpha, Element beta, Double_t &seed)
randomize matrix element values but keep matrix symmetric positive definite
virtual const Int_t * GetColIndexArray() const
TMatrixTSym< Element > & Transpose(const TMatrixTSym< Element > &source)
Transpose a matrix.
Element * New_m(Int_t size)
return data pointer .
TMatrixTSym< Element > & Use(Int_t row_lwb, Int_t row_upb, Element *data)
TMatrixTSym< Element > & Rank1Update(const TVectorT< Element > &v, Element alpha=1.0)
Perform a rank 1 operation on the matrix: A += alpha * v * v^T.
void Minus(const TMatrixTSym< Element > &a, const TMatrixTSym< Element > &b)
Symmetric matrix summation. Create a matrix C such that C = A + B.
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.
TMatrixTSym< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
void TMult(const TMatrixT< Element > &a)
Create a matrix C such that C = A' * A.
void Delete_m(Int_t size, Element *&)
delete data pointer m, if it was assigned on the heap
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 > & 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 TMatrixTBase< Element > & SetMatrixArray(const Element *data, Option_t *option="")
Copy array data to matrix .
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 .
const TMatrixT< Element > EigenVectors(TVectorT< Element > &eigenValues) const
Return a matrix containing the eigen-vectors ordered by descending eigen-values.
TMatrixTSym< Element > & operator-=(Element val)
Subtract val from every element of the matrix.
void Plus(const TMatrixTSym< Element > &a, const TMatrixTSym< Element > &b)
Symmetric matrix summation. Create a matrix C such that C = A + B.
virtual const Int_t * GetRowIndexArray() const
TMatrixTBase< Element > & Apply(const TElementActionT< Element > &action)
Apply action to each matrix element.
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...
virtual TMatrixTBase< Element > & Randomize(Element alpha, Element beta, Double_t &seed)
randomize matrix element values but keep matrix symmetric
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...
TMatrixTSym< Element > & operator=(const TMatrixTSym< Element > &source)
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....
Element * fElements
data container
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...
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])
virtual const Element * GetMatrixArray() const
void TMult(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
Create a matrix C such that C = A' * B.
TObject & operator=(const TObject &rhs)
TObject assignment operator.
void ToUpper()
Change string to upper case.
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
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.
double beta(double x, double y)
Calculates the beta function.
EvaluateInfo init(std::vector< RooRealProxy > parameters, std::vector< ArrayWrapper * > wrappers, std::vector< double * > arrays, size_t begin, size_t batchSize)
RPY_EXPORTED TCppObject_t Allocate(TCppType_t type)
void AMultB(int n, int m, int k, const double *A, const double *B, double *C)
int Invert(LASymMatrix &)
static constexpr double m2
Short_t Max(Short_t a, Short_t b)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Short_t Min(Short_t a, Short_t b)