40template<
class Element>
43 Allocate(no_rows,no_rows,0,0,1);
48template<
class Element>
51 const Int_t no_rows = row_upb-row_lwb+1;
52 Allocate(no_rows,no_rows,row_lwb,row_lwb,1);
64template<
class Element>
67 Allocate(no_rows,no_rows);
68 SetMatrixArray(elements,option);
69 if (!this->IsSymmetric()) {
70 Error(
"TMatrixTSym(Int_t,Element*,Option_t*)",
"matrix not symmetric");
77template<
class Element>
80 const Int_t no_rows = row_upb-row_lwb+1;
81 Allocate(no_rows,no_rows,row_lwb,row_lwb);
82 SetMatrixArray(elements,option);
83 if (!this->IsSymmetric()) {
84 Error(
"TMatrixTSym(Int_t,Int_t,Element*,Option_t*)",
"matrix not symmetric");
90template<
class Element>
103template<
class Element>
123 Transpose(prototype);
133 const Element oldTol = this->SetTol(std::numeric_limits<Element>::min());
135 this->SetTol(oldTol);
145 Error(
"TMatrixTSym(EMatrixCreatorOp1,const TMatrixTSym)",
146 "operation %d not yet implemented", op);
152template<
class Element>
164 Error(
"TMatrixTSym(EMatrixCreatorOp1,const TMatrixT)",
165 "operation %d not yet implemented", op);
171template<
class Element>
180 Allocate(
a.GetNcols(),
a.GetNcols(),
a.GetColLwb(),
a.GetColLwb(),1);
187 Allocate(
a.GetNcols(),
a.GetNcols(),
a.GetColLwb(),
a.GetColLwb(),1);
193 Error(
"TMatrixTSym(EMatrixCreatorOp2)",
"operation %d not yet implemented", op);
199template<
class Element>
205 lazy_constructor.
FillIn(*
this);
206 if (!this->IsSymmetric()) {
207 Error(
"TMatrixTSym(TMatrixTSymLazy)",
"matrix not symmetric");
214template<
class Element>
218 if (
size > this->kSizeMax)
228template<
class Element>
231 if (
size == 0)
return 0;
233 if ( size <= this->kSizeMax )
236 Element *heap =
new Element[
size];
246template<
class Element>
250 if (copySize == 0 || oldp == newp)
253 if ( newSize <= this->kSizeMax && oldSize <= this->kSizeMax ) {
256 for (
Int_t i = copySize-1; i >= 0; i--)
259 for (
Int_t i = 0; i < copySize; i++)
264 memcpy(newp,oldp,copySize*
sizeof(Element));
273template<
class Element>
277 this->fIsOwner =
kTRUE;
278 this->fTol = std::numeric_limits<Element>::epsilon();
286 if (no_rows < 0 || no_cols < 0)
288 Error(
"Allocate",
"no_rows=%d no_cols=%d",no_rows,no_cols);
294 this->fNrows = no_rows;
295 this->fNcols = no_cols;
296 this->fRowLwb = row_lwb;
297 this->fColLwb = col_lwb;
298 this->fNelems = this->fNrows*this->fNcols;
300 if (this->fNelems > 0) {
301 fElements = New_m(this->fNelems);
303 memset(fElements,0,this->fNelems*
sizeof(Element));
311template<
class Element>
316 Error(
"Plus",
"matrices not compatible");
320 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
321 Error(
"Plus",
"this->GetMatrixArray() == a.GetMatrixArray()");
325 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
326 Error(
"Plus",
"this->GetMatrixArray() == b.GetMatrixArray()");
331 const Element * ap =
a.GetMatrixArray();
332 const Element * bp =
b.GetMatrixArray();
333 Element * cp = this->GetMatrixArray();
334 const Element *
const cp_last = cp+this->fNelems;
336 while (cp < cp_last) {
345template<
class Element>
350 Error(
"Minus",
"matrices not compatible");
354 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
355 Error(
"Minus",
"this->GetMatrixArray() == a.GetMatrixArray()");
359 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
360 Error(
"Minus",
"this->GetMatrixArray() == b.GetMatrixArray()");
365 const Element * ap =
a.GetMatrixArray();
366 const Element * bp =
b.GetMatrixArray();
367 Element * cp = this->GetMatrixArray();
368 const Element *
const cp_last = cp+this->fNelems;
370 while (cp < cp_last) {
380template<
class Element>
386 const Element *ap =
a.GetMatrixArray();
387 Element *cp = this->GetMatrixArray();
388 if (
typeid(Element) ==
typeid(
Double_t))
389 cblas_dgemm (CblasRowMajor,CblasTrans,CblasNoTrans,this->fNrows,this->fNcols,
a.GetNrows(),
390 1.0,ap,
a.GetNcols(),ap,
a.GetNcols(),1.0,cp,this->fNcols);
391 else if (
typeid(Element) !=
typeid(
Float_t))
392 cblas_sgemm (CblasRowMajor,CblasTrans,CblasNoTrans,fNrows,fNcols,
a.GetNrows(),
393 1.0,ap,
a.GetNcols(),ap,
a.GetNcols(),1.0,cp,fNcols);
395 Error(
"TMult",
"type %s not implemented in BLAS library",
typeid(Element));
397 const Int_t nb =
a.GetNoElements();
398 const Int_t ncolsa =
a.GetNcols();
399 const Int_t ncolsb = ncolsa;
400 const Element *
const ap =
a.GetMatrixArray();
401 const Element *
const bp = ap;
402 Element * cp = this->GetMatrixArray();
404 const Element *acp0 = ap;
405 while (acp0 < ap+
a.GetNcols()) {
406 for (
const Element *bcp = bp; bcp < bp+ncolsb; ) {
407 const Element *acp = acp0;
409 while (bcp < bp+nb) {
420 R__ASSERT(cp == this->GetMatrixArray()+this->fNelems && acp0 == ap+ncolsa);
428template<
class Element>
434 const Element *ap =
a.GetMatrixArray();
435 Element *cp = this->GetMatrixArray();
436 if (
typeid(Element) ==
typeid(
Double_t))
437 cblas_dsymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
438 ap,
a.GetNcols(),ap,
a.GetNcols(),0.0,cp,this->fNcols);
439 else if (
typeid(Element) !=
typeid(
Float_t))
440 cblas_ssymm (CblasRowMajor,CblasLeft,CblasUpper,fNrows,fNcols,1.0,
441 ap1,
a.GetNcols(),ap,
a.GetNcols(),0.0,cp,fNcols);
443 Error(
"TMult",
"type %s not implemented in BLAS library",
typeid(Element));
445 const Int_t nb =
a.GetNoElements();
446 const Int_t ncolsa =
a.GetNcols();
447 const Int_t ncolsb = ncolsa;
448 const Element *
const ap =
a.GetMatrixArray();
449 const Element *
const bp = ap;
450 Element * cp = this->GetMatrixArray();
452 const Element *acp0 = ap;
453 while (acp0 < ap+
a.GetNcols()) {
454 for (
const Element *bcp = bp; bcp < bp+ncolsb; ) {
455 const Element *acp = acp0;
457 while (bcp < bp+nb) {
468 R__ASSERT(cp == this->GetMatrixArray()+this->fNelems && acp0 == ap+ncolsa);
474template<
class Element>
479 Error(
"Use",
"row_upb=%d < row_lwb=%d",row_upb,row_lwb);
484 this->fNrows = row_upb-row_lwb+1;
485 this->fNcols = this->fNrows;
486 this->fRowLwb = row_lwb;
487 this->fColLwb = row_lwb;
488 this->fNelems = this->fNrows*this->fNcols;
502template<
class Element>
508 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
509 Error(
"GetSub",
"row_lwb out of bounds");
512 if (row_upb < this->fRowLwb || row_upb > this->fRowLwb+this->fNrows-1) {
513 Error(
"GetSub",
"row_upb out of bounds");
516 if (row_upb < row_lwb) {
517 Error(
"GetSub",
"row_upb < row_lwb");
530 row_upb_sub = row_upb-row_lwb;
532 row_lwb_sub = row_lwb;
533 row_upb_sub = row_upb;
536 target.
ResizeTo(row_lwb_sub,row_upb_sub,row_lwb_sub,row_upb_sub);
537 const Int_t nrows_sub = row_upb_sub-row_lwb_sub+1;
540 for (
Int_t irow = 0; irow < nrows_sub; irow++) {
541 for (
Int_t icol = 0; icol < nrows_sub; icol++) {
542 target(irow+row_lwb_sub,icol+row_lwb_sub) = (*this)(row_lwb+irow,row_lwb+icol);
546 const Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNrows+(row_lwb-this->fRowLwb);
549 for (
Int_t irow = 0; irow < nrows_sub; irow++) {
550 const Element *ap_sub = ap;
551 for (
Int_t icol = 0; icol < nrows_sub; icol++) {
568template<
class Element>
574 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
575 Error(
"GetSub",
"row_lwb out of bounds");
578 if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
579 Error(
"GetSub",
"col_lwb out of bounds");
582 if (row_upb < this->fRowLwb || row_upb > this->fRowLwb+this->fNrows-1) {
583 Error(
"GetSub",
"row_upb out of bounds");
586 if (col_upb < this->fColLwb || col_upb > this->fColLwb+this->fNcols-1) {
587 Error(
"GetSub",
"col_upb out of bounds");
590 if (row_upb < row_lwb || col_upb < col_lwb) {
591 Error(
"GetSub",
"row_upb < row_lwb || col_upb < col_lwb");
600 const Int_t row_lwb_sub = (shift) ? 0 : row_lwb;
601 const Int_t row_upb_sub = (shift) ? row_upb-row_lwb : row_upb;
602 const Int_t col_lwb_sub = (shift) ? 0 : col_lwb;
603 const Int_t col_upb_sub = (shift) ? col_upb-col_lwb : col_upb;
605 target.
ResizeTo(row_lwb_sub,row_upb_sub,col_lwb_sub,col_upb_sub);
606 const Int_t nrows_sub = row_upb_sub-row_lwb_sub+1;
607 const Int_t ncols_sub = col_upb_sub-col_lwb_sub+1;
610 for (
Int_t irow = 0; irow < nrows_sub; irow++) {
611 for (
Int_t icol = 0; icol < ncols_sub; icol++) {
612 target(irow+row_lwb_sub,icol+col_lwb_sub) = (*this)(row_lwb+irow,col_lwb+icol);
616 const Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNcols+(col_lwb-this->fColLwb);
619 for (
Int_t irow = 0; irow < nrows_sub; irow++) {
620 const Element *ap_sub = ap;
621 for (
Int_t icol = 0; icol < ncols_sub; icol++) {
635template<
class Element>
643 Error(
"SetSub",
"source matrix is not symmetric");
646 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
647 Error(
"SetSub",
"row_lwb outof bounds");
650 if (row_lwb+source.
GetNrows() > this->fRowLwb+this->fNrows) {
651 Error(
"SetSub",
"source matrix too large");
660 for (
Int_t irow = 0; irow < nRows_source; irow++) {
661 for (
Int_t icol = 0; icol < nRows_source; icol++) {
662 (*this)(row_lwb+irow,row_lwb+icol) = source(rowlwb_s+irow,rowlwb_s+icol);
667 Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNrows+(row_lwb-this->fRowLwb);
669 for (
Int_t irow = 0; irow < nRows_source; irow++) {
670 Element *ap_sub = ap;
671 for (
Int_t icol = 0; icol < nRows_source; icol++) {
685template<
class Element>
692 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
693 Error(
"SetSub",
"row_lwb out of bounds");
696 if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
697 Error(
"SetSub",
"col_lwb out of bounds");
701 if (row_lwb+source.
GetNrows() > this->fRowLwb+this->fNrows || col_lwb+source.
GetNcols() > this->fRowLwb+this->fNrows) {
702 Error(
"SetSub",
"source matrix too large");
705 if (col_lwb+source.
GetNcols() > this->fRowLwb+this->fNrows || row_lwb+source.
GetNrows() > this->fRowLwb+this->fNrows) {
706 Error(
"SetSub",
"source matrix too large");
716 if (row_lwb >= col_lwb) {
719 for (irow = 0; irow < nRows_source; irow++) {
720 for (
Int_t icol = 0; col_lwb+icol <= row_lwb+irow &&
721 icol < nCols_source; icol++) {
722 (*this)(row_lwb+irow,col_lwb+icol) = source(irow+rowlwb_s,icol+collwb_s);
727 for (irow = 0; irow < nCols_source; irow++) {
728 for (
Int_t icol = nRows_source-1; row_lwb+icol > irow+col_lwb &&
730 (*this)(col_lwb+irow,row_lwb+icol) = source(icol+rowlwb_s,irow+collwb_s);
742template<
class Element>
746 if (!this->IsSymmetric()) {
747 Error(
"SetMatrixArray",
"Matrix is not symmetric after Set");
755template<
class Element>
758 if (row_shift != col_shift) {
759 Error(
"Shift",
"row_shift != col_shift");
770template<
class Element>
774 if (!this->fIsOwner) {
775 Error(
"ResizeTo(Int_t,Int_t)",
"Not owner of data array,cannot resize");
779 if (nrows != ncols) {
780 Error(
"ResizeTo(Int_t,Int_t)",
"nrows != ncols");
784 if (this->fNelems > 0) {
785 if (this->fNrows == nrows && this->fNcols == ncols)
787 else if (nrows == 0 || ncols == 0) {
788 this->fNrows = nrows; this->fNcols = ncols;
793 Element *elements_old = GetMatrixArray();
794 const Int_t nelems_old = this->fNelems;
795 const Int_t nrows_old = this->fNrows;
796 const Int_t ncols_old = this->fNcols;
798 Allocate(nrows,ncols);
801 Element *elements_new = GetMatrixArray();
804 if (this->fNelems > this->kSizeMax || nelems_old > this->kSizeMax)
805 memset(elements_new,0,this->fNelems*
sizeof(Element));
806 else if (this->fNelems > nelems_old)
807 memset(elements_new+nelems_old,0,(this->fNelems-nelems_old)*
sizeof(Element));
813 const Int_t nelems_new = this->fNelems;
814 if (ncols_old < this->fNcols) {
815 for (
Int_t i = nrows_copy-1; i >= 0; i--) {
816 Memcpy_m(elements_new+i*this->fNcols,elements_old+i*ncols_old,ncols_copy,
817 nelems_new,nelems_old);
818 if (this->fNelems <= this->kSizeMax && nelems_old <= this->kSizeMax)
819 memset(elements_new+i*this->fNcols+ncols_copy,0,(this->fNcols-ncols_copy)*
sizeof(Element));
822 for (
Int_t i = 0; i < nrows_copy; i++)
823 Memcpy_m(elements_new+i*this->fNcols,elements_old+i*ncols_old,ncols_copy,
824 nelems_new,nelems_old);
827 Delete_m(nelems_old,elements_old);
829 Allocate(nrows,ncols,0,0,1);
840template<
class Element>
845 if (!this->fIsOwner) {
846 Error(
"ResizeTo(Int_t,Int_t,Int_t,Int_t)",
"Not owner of data array,cannot resize");
850 if (row_lwb != col_lwb) {
851 Error(
"ResizeTo(Int_t,Int_t,Int_t,Int_t)",
"row_lwb != col_lwb");
854 if (row_upb != col_upb) {
855 Error(
"ResizeTo(Int_t,Int_t,Int_t,Int_t)",
"row_upb != col_upb");
859 const Int_t new_nrows = row_upb-row_lwb+1;
860 const Int_t new_ncols = col_upb-col_lwb+1;
862 if (this->fNelems > 0) {
864 if (this->fNrows == new_nrows && this->fNcols == new_ncols &&
865 this->fRowLwb == row_lwb && this->fColLwb == col_lwb)
867 else if (new_nrows == 0 || new_ncols == 0) {
868 this->fNrows = new_nrows; this->fNcols = new_ncols;
869 this->fRowLwb = row_lwb; this->fColLwb = col_lwb;
874 Element *elements_old = GetMatrixArray();
875 const Int_t nelems_old = this->fNelems;
876 const Int_t nrows_old = this->fNrows;
877 const Int_t ncols_old = this->fNcols;
878 const Int_t rowLwb_old = this->fRowLwb;
879 const Int_t colLwb_old = this->fColLwb;
881 Allocate(new_nrows,new_ncols,row_lwb,col_lwb);
884 Element *elements_new = GetMatrixArray();
887 if (this->fNelems > this->kSizeMax || nelems_old > this->kSizeMax)
888 memset(elements_new,0,this->fNelems*
sizeof(Element));
889 else if (this->fNelems > nelems_old)
890 memset(elements_new+nelems_old,0,(this->fNelems-nelems_old)*
sizeof(Element));
895 const Int_t rowUpb_copy =
TMath::Min(this->fRowLwb+this->fNrows-1,rowLwb_old+nrows_old-1);
896 const Int_t colUpb_copy =
TMath::Min(this->fColLwb+this->fNcols-1,colLwb_old+ncols_old-1);
898 const Int_t nrows_copy = rowUpb_copy-rowLwb_copy+1;
899 const Int_t ncols_copy = colUpb_copy-colLwb_copy+1;
901 if (nrows_copy > 0 && ncols_copy > 0) {
902 const Int_t colOldOff = colLwb_copy-colLwb_old;
903 const Int_t colNewOff = colLwb_copy-this->fColLwb;
904 if (ncols_old < this->fNcols) {
905 for (
Int_t i = nrows_copy-1; i >= 0; i--) {
906 const Int_t iRowOld = rowLwb_copy+i-rowLwb_old;
907 const Int_t iRowNew = rowLwb_copy+i-this->fRowLwb;
908 Memcpy_m(elements_new+iRowNew*this->fNcols+colNewOff,
909 elements_old+iRowOld*ncols_old+colOldOff,ncols_copy,this->fNelems,nelems_old);
910 if (this->fNelems <= this->kSizeMax && nelems_old <= this->kSizeMax)
911 memset(elements_new+iRowNew*this->fNcols+colNewOff+ncols_copy,0,
912 (this->fNcols-ncols_copy)*
sizeof(Element));
915 for (
Int_t i = 0; i < nrows_copy; i++) {
916 const Int_t iRowOld = rowLwb_copy+i-rowLwb_old;
917 const Int_t iRowNew = rowLwb_copy+i-this->fRowLwb;
918 Memcpy_m(elements_new+iRowNew*this->fNcols+colNewOff,
919 elements_old+iRowOld*ncols_old+colOldOff,ncols_copy,this->fNelems,nelems_old);
924 Delete_m(nelems_old,elements_old);
926 Allocate(new_nrows,new_ncols,row_lwb,col_lwb,1);
934template<
class Element>
946template<
class Element>
960template<
class Element>
967 Element *p2 = this->GetMatrixArray();
968 for (
Int_t i = 0; i < this->GetNoElements(); i++)
978template<
class Element>
987 Element *pM = this->GetMatrixArray();
989 Error(
"InvertFast",
"matrix is singular");
999 TMatrixTSymCramerInv::Inv2x2<Element>(*
this,det);
1004 TMatrixTSymCramerInv::Inv3x3<Element>(*
this,det);
1009 TMatrixTSymCramerInv::Inv4x4<Element>(*
this,det);
1014 TMatrixTSymCramerInv::Inv5x5<Element>(*
this,det);
1019 TMatrixTSymCramerInv::Inv6x6<Element>(*
this,det);
1028 Element *p2 = this->GetMatrixArray();
1029 for (
Int_t i = 0; i < this->GetNoElements(); i++)
1040template<
class Element>
1049 Error(
"Transpose",
"matrix has wrong shape");
1062template<
class Element>
1068 if (
v.GetNoElements() < this->fNrows) {
1069 Error(
"Rank1Update",
"vector too short");
1074 const Element *
const pv =
v.GetMatrixArray();
1075 Element *trp = this->GetMatrixArray();
1077 for (
Int_t i = 0; i < this->fNrows; i++) {
1079 tcp += i*this->fNcols;
1080 const Element tmp = alpha*pv[i];
1081 for (
Int_t j = i; j < this->fNcols; j++) {
1082 if (j > i) *tcp += tmp*pv[j];
1083 *trp++ += tmp*pv[j];
1084 tcp += this->fNcols;
1086 tcp -= this->fNelems-1;
1098template<
class Element>
1104 if (this->fNcols !=
b.GetNcols() || this->fColLwb !=
b.GetColLwb()) {
1105 Error(
"Similarity(const TMatrixT &)",
"matrices incompatible");
1110 const Int_t ncolsa = this->fNcols;
1111 const Int_t nb =
b.GetNoElements();
1112 const Int_t nrowsb =
b.GetNrows();
1113 const Int_t ncolsb =
b.GetNcols();
1115 const Element *
const bp =
b.GetMatrixArray();
1117 Element work[kWorkMax];
1119 Element *bap = work;
1120 if (nrowsb*ncolsa > kWorkMax) {
1121 isAllocated =
kTRUE;
1122 bap =
new Element[nrowsb*ncolsa];
1125 AMultB(bp,nb,ncolsb,this->fElements,this->fNelems,this->fNcols,bap);
1127 if (nrowsb != this->fNrows)
1128 this->ResizeTo(nrowsb,nrowsb);
1131 Element *cp = this->GetMatrixArray();
1132 if (
typeid(Element) ==
typeid(
Double_t))
1133 cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasTrans,this->fNrows,this->fNcols,ba.GetNcols(),
1134 1.0,bap,ba.GetNcols(),bp,
b.GetNcols(),1.0,cp,this->fNcols);
1135 else if (
typeid(Element) !=
typeid(
Float_t))
1136 cblas_sgemm (CblasRowMajor,CblasNoTrans,CblasTrans,this->fNrows,this->fNcols,ba.GetNcols(),
1137 1.0,bap,ba.GetNcols(),bp,
b.GetNcols(),1.0,cp,this->fNcols);
1139 Error(
"Similarity",
"type %s not implemented in BLAS library",
typeid(Element));
1141 const Int_t nba = nrowsb*ncolsa;
1142 const Int_t ncolsba = ncolsa;
1143 const Element * bi1p = bp;
1144 Element * cp = this->GetMatrixArray();
1145 Element *
const cp0 = cp;
1148 const Element *barp0 = bap;
1149 while (barp0 < bap+nba) {
1150 const Element *brp0 = bi1p;
1151 while (brp0 < bp+nb) {
1152 const Element *barp = barp0;
1153 const Element *brp = brp0;
1155 while (brp < brp0+ncolsb)
1156 cij += *barp++ * *brp++;
1165 R__ASSERT(cp == cp0+this->fNelems+ishift && barp0 == bap+nba);
1168 for (
Int_t irow = 0; irow < this->fNrows; irow++) {
1169 const Int_t rowOff1 = irow*this->fNrows;
1170 for (
Int_t icol = 0; icol < irow; icol++) {
1171 const Int_t rowOff2 = icol*this->fNrows;
1172 cp[rowOff1+icol] = cp[rowOff2+irow];
1189template<
class Element>
1195 if (this->fNcols !=
b.GetNcols() || this->fColLwb !=
b.GetColLwb()) {
1196 Error(
"Similarity(const TMatrixTSym &)",
"matrices incompatible");
1202 const Int_t nrowsb =
b.GetNrows();
1203 const Int_t ncolsa = this->GetNcols();
1205 Element work[kWorkMax];
1207 Element *abtp = work;
1208 if (this->fNcols > kWorkMax) {
1209 isAllocated =
kTRUE;
1210 abtp =
new Element[this->fNcols];
1215 const Element *bp =
b.GetMatrixArray();
1216 Element *cp = this->GetMatrixArray();
1217 if (
typeid(Element) ==
typeid(
Double_t))
1218 cblas_dsymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
1219 bp,
b.GetNcols(),abtp,abt.
GetNcols(),0.0,cp,this->fNcols);
1220 else if (
typeid(Element) !=
typeid(
Float_t))
1221 cblas_ssymm (CblasRowMajor,CblasLeft,CblasUpper,this->fNrows,this->fNcols,1.0,
1222 bp,
b.GetNcols(),abtp,abt.
GetNcols(),0.0,cp,this->fNcols);
1224 Error(
"Similarity",
"type %s not implemented in BLAS library",
typeid(Element));
1229 const Int_t ncolsa = this->GetNcols();
1230 const Int_t nb =
b.GetNoElements();
1231 const Int_t nrowsb =
b.GetNrows();
1232 const Int_t ncolsb =
b.GetNcols();
1234 const Element *
const bp =
b.GetMatrixArray();
1236 Element work[kWorkMax];
1238 Element *bap = work;
1239 if (nrowsb*ncolsa > kWorkMax) {
1240 isAllocated =
kTRUE;
1241 bap =
new Element[nrowsb*ncolsa];
1244 AMultB(bp,nb,ncolsb,this->fElements,this->fNelems,this->fNcols,bap);
1246 const Int_t nba = nrowsb*ncolsa;
1247 const Int_t ncolsba = ncolsa;
1248 const Element * bi1p = bp;
1249 Element * cp = this->GetMatrixArray();
1250 Element *
const cp0 = cp;
1253 const Element *barp0 = bap;
1254 while (barp0 < bap+nba) {
1255 const Element *brp0 = bi1p;
1256 while (brp0 < bp+nb) {
1257 const Element *barp = barp0;
1258 const Element *brp = brp0;
1260 while (brp < brp0+ncolsb)
1261 cij += *barp++ * *brp++;
1270 R__ASSERT(cp == cp0+this->fNelems+ishift && barp0 == bap+nba);
1273 for (
Int_t irow = 0; irow < this->fNrows; irow++) {
1274 const Int_t rowOff1 = irow*this->fNrows;
1275 for (
Int_t icol = 0; icol < irow; icol++) {
1276 const Int_t rowOff2 = icol*this->fNrows;
1277 cp[rowOff1+icol] = cp[rowOff2+irow];
1291template<
class Element>
1297 if (this->fNcols !=
v.GetNrows() || this->fColLwb !=
v.GetLwb()) {
1298 Error(
"Similarity(const TVectorT &)",
"vector and matrix incompatible");
1303 const Element *mp = this->GetMatrixArray();
1304 const Element *vp =
v.GetMatrixArray();
1307 const Element *
const vp_first = vp;
1308 const Element *
const vp_last = vp+
v.GetNrows();
1309 while (vp < vp_last) {
1311 for (
const Element *sp = vp_first; sp < vp_last; )
1312 sum2 += *mp++ * *sp++;
1313 sum1 += sum2 * *vp++;
1316 R__ASSERT(mp == this->GetMatrixArray()+this->GetNoElements());
1326template<
class Element>
1332 if (this->fNrows !=
b.GetNrows() || this->fRowLwb !=
b.GetRowLwb()) {
1333 Error(
"SimilarityT(const TMatrixT &)",
"matrices incompatible");
1338 const Int_t ncolsb =
b.GetNcols();
1339 const Int_t ncolsa = this->GetNcols();
1341 Element work[kWorkMax];
1343 Element *btap = work;
1344 if (ncolsb*ncolsa > kWorkMax) {
1345 isAllocated =
kTRUE;
1346 btap =
new Element[ncolsb*ncolsa];
1352 if (ncolsb != this->fNcols)
1353 this->ResizeTo(ncolsb,ncolsb);
1356 const Element *bp =
b.GetMatrixArray();
1357 Element *cp = this->GetMatrixArray();
1358 if (
typeid(Element) ==
typeid(
Double_t))
1359 cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasNoTrans,this->fNrows,this->fNcols,bta.
GetNcols(),
1360 1.0,btap,bta.
GetNcols(),bp,
b.GetNcols(),1.0,cp,this->fNcols);
1361 else if (
typeid(Element) !=
typeid(
Float_t))
1362 cblas_sgemm (CblasRowMajor,CblasNoTrans,CblasNoTrans,this->fNrows,this->fNcols,bta.
GetNcols(),
1363 1.0,btap,bta.
GetNcols(),bp,
b.GetNcols(),1.0,cp,this->fNcols);
1365 Error(
"similarityT",
"type %s not implemented in BLAS library",
typeid(Element));
1368 const Int_t nb =
b.GetNoElements();
1370 const Element *
const bp =
b.GetMatrixArray();
1371 Element * cp = this->GetMatrixArray();
1372 Element *
const cp0 = cp;
1375 const Element *btarp0 = btap;
1376 const Element *bcp0 = bp;
1377 while (btarp0 < btap+nbta) {
1378 for (
const Element *bcp = bcp0; bcp < bp+ncolsb; ) {
1379 const Element *btarp = btarp0;
1381 while (bcp < bp+nb) {
1382 cij += *btarp++ * *bcp;
1393 R__ASSERT(cp == cp0+this->fNelems+ishift && btarp0 == btap+nbta);
1396 for (
Int_t irow = 0; irow < this->fNrows; irow++) {
1397 const Int_t rowOff1 = irow*this->fNrows;
1398 for (
Int_t icol = 0; icol < irow; icol++) {
1399 const Int_t rowOff2 = icol*this->fNrows;
1400 cp[rowOff1+icol] = cp[rowOff2+irow];
1413template<
class Element>
1417 Error(
"operator=",
"matrices not compatible");
1423 memcpy(this->GetMatrixArray(),source.
fElements,this->fNelems*
sizeof(Element));
1430template<
class Element>
1435 if (lazy_constructor.
fRowUpb != this->GetRowUpb() ||
1436 lazy_constructor.
fRowLwb != this->GetRowLwb()) {
1437 Error(
"operator=(const TMatrixTSymLazy&)",
"matrix is incompatible with "
1438 "the assigned Lazy matrix");
1442 lazy_constructor.
FillIn(*
this);
1449template<
class Element>
1454 Element *ep = fElements;
1455 const Element *
const ep_last = ep+this->fNelems;
1456 while (ep < ep_last)
1465template<
class Element>
1470 Element *ep = fElements;
1471 const Element *
const ep_last = ep+this->fNelems;
1472 while (ep < ep_last)
1481template<
class Element>
1486 Element *ep = fElements;
1487 const Element *
const ep_last = ep+this->fNelems;
1488 while (ep < ep_last)
1497template<
class Element>
1502 Element *ep = fElements;
1503 const Element *
const ep_last = ep+this->fNelems;
1504 while (ep < ep_last)
1513template<
class Element>
1517 Error(
"operator+=",
"matrices not compatible");
1522 Element *tp = this->GetMatrixArray();
1523 const Element *
const tp_last = tp+this->fNelems;
1524 while (tp < tp_last)
1533template<
class Element>
1537 Error(
"operator-=",
"matrices not compatible");
1542 Element *tp = this->GetMatrixArray();
1543 const Element *
const tp_last = tp+this->fNelems;
1544 while (tp < tp_last)
1552template<
class Element>
1558 Element *trp = this->GetMatrixArray();
1560 for (
Int_t i = 0; i < this->fNrows; i++) {
1562 tcp += i*this->fNcols;
1563 for (
Int_t j = i; j < this->fNcols; j++) {
1565 if (j > i) *tcp = val;
1567 tcp += this->fNcols;
1569 tcp -= this->fNelems-1;
1579template<
class Element>
1585 Element *trp = this->GetMatrixArray();
1587 for (
Int_t i = 0; i < this->fNrows; i++) {
1588 action.
fI = i+this->fRowLwb;
1590 tcp += i*this->fNcols;
1591 for (
Int_t j = i; j < this->fNcols; j++) {
1592 action.
fJ = j+this->fColLwb;
1594 if (j > i) *tcp = val;
1596 tcp += this->fNcols;
1598 tcp -= this->fNelems-1;
1607template<
class Element>
1612 if (this->fNrows != this->fNcols || this->fRowLwb != this->fColLwb) {
1613 Error(
"Randomize(Element,Element,Element &",
"matrix should be square");
1618 const Element scale = beta-alpha;
1619 const Element shift = alpha/scale;
1621 Element *ep = GetMatrixArray();
1622 for (
Int_t i = 0; i < this->fNrows; i++) {
1623 const Int_t off = i*this->fNcols;
1624 for (
Int_t j = 0; j <= i; j++) {
1625 ep[off+j] = scale*(
Drand(seed)+shift);
1627 ep[j*this->fNcols+i] = ep[off+j];
1638template<
class Element>
1643 if (this->fNrows != this->fNcols || this->fRowLwb != this->fColLwb) {
1644 Error(
"RandomizeSym(Element,Element,Element &",
"matrix should be square");
1649 const Element scale = beta-alpha;
1650 const Element shift = alpha/scale;
1652 Element *ep = GetMatrixArray();
1654 for (i = 0; i < this->fNrows; i++) {
1655 const Int_t off = i*this->fNcols;
1656 for (
Int_t j = 0; j <= i; j++)
1657 ep[off+j] = scale*(
Drand(seed)+shift);
1660 for (i = this->fNrows-1; i >= 0; i--) {
1661 const Int_t off1 = i*this->fNcols;
1662 for (
Int_t j = i; j >= 0; j--) {
1663 const Int_t off2 = j*this->fNcols;
1664 ep[off1+j] *= ep[off2+j];
1665 for (
Int_t k = j-1; k >= 0; k--) {
1666 ep[off1+j] += ep[off1+k]*ep[off2+k];
1669 ep[off2+i] = ep[off1+j];
1680template<
class Element>
1685 eigenValues.
ResizeTo(this->fNrows);
1693template<
class Element>
1703template<
class Element>
1713template<
class Element>
1723template<
class Element>
1731template<
class Element>
1741template<
class Element>
1751template<
class Element>
1754 return Element(-1.0)*
operator-(source1,val);
1759template<
class Element>
1769template<
class Element>
1778template<
class Element>
1784 Error(
"operator&&(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1794 while (tp < tp_last)
1795 *tp++ = (*sp1++ != 0.0 && *sp2++ != 0.0);
1803template<
class Element>
1809 Error(
"operator||(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1819 while (tp < tp_last)
1820 *tp++ = (*sp1++ != 0.0 || *sp2++ != 0.0);
1828template<
class Element>
1834 Error(
"operator>(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1844 while (tp < tp_last) {
1845 *tp++ = (*sp1) > (*sp2); sp1++; sp2++;
1854template<
class Element>
1860 Error(
"operator>=(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1870 while (tp < tp_last) {
1871 *tp++ = (*sp1) >= (*sp2); sp1++; sp2++;
1880template<
class Element>
1886 Error(
"operator<=(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1896 while (tp < tp_last) {
1897 *tp++ = (*sp1) <= (*sp2); sp1++; sp2++;
1906template<
class Element>
1912 Error(
"operator<(const TMatrixTSym&,const TMatrixTSym&)",
"matrices not compatible");
1922 while (tp < tp_last) {
1923 *tp++ = (*sp1) < (*sp2); sp1++; sp2++;
1932template<
class Element>
1936 ::Error(
"Add",
"matrices not compatible");
1946 for (
Int_t i = 0; i < nrows; i++) {
1950 for (
Int_t j = i; j < ncols; j++) {
1951 const Element tmp = scalar * *sp++;
1952 if (j > i) *tcp += tmp;
1965template<
class Element>
1969 ::Error(
"ElementMult",
"matrices not compatible");
1979 for (
Int_t i = 0; i < nrows; i++) {
1983 for (
Int_t j = i; j < ncols; j++) {
1984 if (j > i) *tcp *= *sp;
1997template<
class Element>
2001 ::Error(
"ElementDiv",
"matrices not compatible");
2011 for (
Int_t i = 0; i < nrows; i++) {
2015 for (
Int_t j = i; j < ncols; j++) {
2017 if (j > i) *tcp /= *sp;
2022 Error(
"ElementDiv",
"source (%d,%d) is zero",irow,icol);
2036template<
class Element>
2044 fElements =
new Element[this->fNelems];
2046 for (i = 0; i < this->fNrows; i++) {
2047 R__b.
ReadFastArray(fElements+i*this->fNcols+i,this->fNcols-i);
2050 for (i = 0; i < this->fNrows; i++) {
2051 for (
Int_t j = 0; j < i; j++) {
2052 fElements[i*this->fNcols+j] = fElements[j*this->fNrows+i];
2055 if (this->fNelems <= this->kSizeMax) {
2056 memcpy(fDataStack,fElements,this->fNelems*
sizeof(Element));
2057 delete [] fElements;
2058 fElements = fDataStack;
2063 for (
Int_t i = 0; i < this->fNrows; i++) {
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define templateClassImp(name)
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
R__EXTERN Int_t gMatrixCheck
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)
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 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.
Bool_t operator==(const TMatrixTSym< Element > &m1, const TMatrixTSym< Element > &m2)
Check to see if two matrices are identical.
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] .
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)
void AMultB(const Element *const ap, Int_t na, Int_t ncolsa, const Element *const bp, Int_t nb, Int_t ncolsb, Element *cp)
Elementary routine to calculate matrix multiplication A*B.
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 .