39template<
class Element>
42 Allocate(nrows,ncols,0,0,1);
48template<
class Element>
51 Allocate(row_upb-row_lwb+1,col_upb-col_lwb+1,row_lwb,col_lwb,1);
62template<
class Element>
65 Allocate(no_rows,no_cols);
72template<
class Element>
74 const Element *elements,
Option_t *option)
76 Allocate(row_upb-row_lwb+1,col_upb-col_lwb+1,row_lwb,col_lwb);
83template<
class Element>
94template<
class Element>
105template<
class Element>
118template<
class Element>
138 Transpose(prototype);
148 const Element oldTol = this->SetTol(std::numeric_limits<Element>::min());
150 this->SetTol(oldTol);
156 TMult(prototype,prototype);
160 Error(
"TMatrixT(EMatrixCreatorOp1)",
"operation %d not yet implemented", op);
169template<
class Element>
177 Allocate(
a.GetNrows(),
b.GetNcols(),
a.GetRowLwb(),
b.GetColLwb(),1);
182 Allocate(
a.GetNcols(),
b.GetNcols(),
a.GetColLwb(),
b.GetColLwb(),1);
187 Allocate(
a.GetNrows(),
b.GetNrows(),
a.GetRowLwb(),
b.GetRowLwb(),1);
193 Allocate(
a.GetNrows(),
a.GetNcols(),
a.GetRowLwb(),
a.GetColLwb(),1);
195 const Element oldTol = this->SetTol(std::numeric_limits<Element>::min());
197 this->SetTol(oldTol);
204 Allocate(
a.GetNrows(),
a.GetNcols(),
a.GetRowLwb(),
a.GetColLwb(),1);
211 Allocate(
a.GetNrows(),
a.GetNcols(),
a.GetRowLwb(),
a.GetColLwb(),1);
217 Error(
"TMatrixT(EMatrixCreatorOp2)",
"operation %d not yet implemented", op);
226template<
class Element>
234 Allocate(
a.GetNrows(),
b.GetNcols(),
a.GetRowLwb(),
b.GetColLwb(),1);
239 Allocate(
a.GetNcols(),
b.GetNcols(),
a.GetColLwb(),
b.GetColLwb(),1);
244 Allocate(
a.GetNrows(),
b.GetNrows(),
a.GetRowLwb(),
b.GetRowLwb(),1);
250 Allocate(
a.GetNrows(),
a.GetNcols(),
a.GetRowLwb(),
a.GetColLwb(),1);
252 const Element oldTol = this->SetTol(std::numeric_limits<Element>::min());
254 this->SetTol(oldTol);
261 Allocate(
a.GetNrows(),
a.GetNcols(),
a.GetRowLwb(),
a.GetColLwb(),1);
268 Allocate(
a.GetNrows(),
a.GetNcols(),
a.GetRowLwb(),
a.GetColLwb(),1);
274 Error(
"TMatrixT(EMatrixCreatorOp2)",
"operation %d not yet implemented", op);
283template<
class Element>
291 Allocate(
a.GetNrows(),
b.GetNcols(),
a.GetRowLwb(),
b.GetColLwb(),1);
296 Allocate(
a.GetNcols(),
b.GetNcols(),
a.GetColLwb(),
b.GetColLwb(),1);
301 Allocate(
a.GetNrows(),
b.GetNrows(),
a.GetRowLwb(),
b.GetRowLwb(),1);
307 Allocate(
a.GetNrows(),
a.GetNcols(),
a.GetRowLwb(),
a.GetColLwb(),1);
309 const Element oldTol = this->SetTol(std::numeric_limits<Element>::min());
311 this->SetTol(oldTol);
318 Allocate(
a.GetNrows(),
a.GetNcols(),
a.GetRowLwb(),
a.GetColLwb(),1);
325 Allocate(
a.GetNrows(),
a.GetNcols(),
a.GetRowLwb(),
a.GetColLwb(),1);
331 Error(
"TMatrixT(EMatrixCreatorOp2)",
"operation %d not yet implemented", op);
340template<
class Element>
348 Allocate(
a.GetNrows(),
b.GetNcols(),
a.GetRowLwb(),
b.GetColLwb(),1);
353 Allocate(
a.GetNcols(),
b.GetNcols(),
a.GetColLwb(),
b.GetColLwb(),1);
358 Allocate(
a.GetNrows(),
b.GetNrows(),
a.GetRowLwb(),
b.GetRowLwb(),1);
364 Allocate(
a.GetNrows(),
a.GetNcols(),
a.GetRowLwb(),
a.GetColLwb(),1);
366 const Element oldTol = this->SetTol(std::numeric_limits<Element>::min());
368 this->SetTol(oldTol);
375 Allocate(
a.GetNrows(),
a.GetNcols(),
a.GetRowLwb(),
a.GetColLwb(),1);
382 Allocate(
a.GetNrows(),
a.GetNcols(),
a.GetRowLwb(),
a.GetColLwb(),1);
388 Error(
"TMatrixT(EMatrixCreatorOp2)",
"operation %d not yet implemented", op);
395template<
class Element>
401 lazy_constructor.
FillIn(*
this);
407template<
class Element>
411 if (
size > this->kSizeMax)
421template<
class Element>
424 if (
size == 0)
return 0;
426 if ( size <= this->kSizeMax )
429 Element *heap =
new Element[
size];
439template<
class Element>
443 if (copySize == 0 || oldp == newp)
446 if ( newSize <= this->kSizeMax && oldSize <= this->kSizeMax ) {
449 for (
Int_t i = copySize-1; i >= 0; i--)
452 for (
Int_t i = 0; i < copySize; i++)
457 memcpy(newp,oldp,copySize*
sizeof(Element));
466template<
class Element>
470 this->fIsOwner =
kTRUE;
471 this->fTol = std::numeric_limits<Element>::epsilon();
479 if (no_rows < 0 || no_cols < 0)
481 Error(
"Allocate",
"no_rows=%d no_cols=%d",no_rows,no_cols);
487 this->fNrows = no_rows;
488 this->fNcols = no_cols;
489 this->fRowLwb = row_lwb;
490 this->fColLwb = col_lwb;
491 this->fNelems = this->fNrows*this->fNcols;
494 if( ((
Long64_t)this->fNrows)*this->fNcols != this->fNelems )
496 Error(
"Allocate",
"too large: no_rows=%d no_cols=%d",no_rows,no_cols);
501 if (this->fNelems > 0) {
502 fElements = New_m(this->fNelems);
504 memset(fElements,0,this->fNelems*
sizeof(Element));
512template<
class Element>
517 Error(
"Plus",
"matrices not compatible");
521 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
522 Error(
"Plus",
"this->GetMatrixArray() == a.GetMatrixArray()");
526 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
527 Error(
"Plus",
"this->GetMatrixArray() == b.GetMatrixArray()");
532 const Element * ap =
a.GetMatrixArray();
533 const Element * bp =
b.GetMatrixArray();
534 Element * cp = this->GetMatrixArray();
535 const Element *
const cp_last = cp+this->fNelems;
537 while (cp < cp_last) {
546template<
class Element>
551 Error(
"Plus",
"matrices not compatible");
555 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
556 Error(
"Plus",
"this->GetMatrixArray() == a.GetMatrixArray()");
560 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
561 Error(
"Plus",
"this->GetMatrixArray() == b.GetMatrixArray()");
566 const Element * ap =
a.GetMatrixArray();
567 const Element * bp =
b.GetMatrixArray();
568 Element * cp = this->GetMatrixArray();
569 const Element *
const cp_last = cp+this->fNelems;
571 while (cp < cp_last) {
580template<
class Element>
585 Error(
"Minus",
"matrices not compatible");
589 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
590 Error(
"Minus",
"this->GetMatrixArray() == a.GetMatrixArray()");
594 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
595 Error(
"Minus",
"this->GetMatrixArray() == b.GetMatrixArray()");
600 const Element * ap =
a.GetMatrixArray();
601 const Element * bp =
b.GetMatrixArray();
602 Element * cp = this->GetMatrixArray();
603 const Element *
const cp_last = cp+this->fNelems;
605 while (cp < cp_last) {
614template<
class Element>
619 Error(
"Minus",
"matrices not compatible");
623 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
624 Error(
"Minus",
"this->GetMatrixArray() == a.GetMatrixArray()");
628 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
629 Error(
"Minus",
"this->GetMatrixArray() == b.GetMatrixArray()");
634 const Element * ap =
a.GetMatrixArray();
635 const Element * bp =
b.GetMatrixArray();
636 Element * cp = this->GetMatrixArray();
637 const Element *
const cp_last = cp+this->fNelems;
639 while (cp < cp_last) {
648template<
class Element>
652 if (
a.GetNcols() !=
b.GetNrows() ||
a.GetColLwb() !=
b.GetRowLwb()) {
653 Error(
"Mult",
"A rows and B columns incompatible");
657 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
658 Error(
"Mult",
"this->GetMatrixArray() == a.GetMatrixArray()");
662 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
663 Error(
"Mult",
"this->GetMatrixArray() == b.GetMatrixArray()");
669 const Element *ap =
a.GetMatrixArray();
670 const Element *bp =
b.GetMatrixArray();
671 Element *cp = this->GetMatrixArray();
672 if (
typeid(Element) ==
typeid(
Double_t))
673 cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasNoTrans,fNrows,fNcols,
a.GetNcols(),
674 1.0,ap,
a.GetNcols(),bp,
b.GetNcols(),1.0,cp,fNcols);
675 else if (
typeid(Element) !=
typeid(
Float_t))
676 cblas_sgemm (CblasRowMajor,CblasNoTrans,CblasNoTrans,fNrows,fNcols,
a.GetNcols(),
677 1.0,ap,
a.GetNcols(),bp,
b.GetNcols(),1.0,cp,fNcols);
679 Error(
"Mult",
"type %s not implemented in BLAS library",
typeid(Element));
681 const Int_t na =
a.GetNoElements();
682 const Int_t nb =
b.GetNoElements();
683 const Int_t ncolsa =
a.GetNcols();
684 const Int_t ncolsb =
b.GetNcols();
685 const Element *
const ap =
a.GetMatrixArray();
686 const Element *
const bp =
b.GetMatrixArray();
687 Element * cp = this->GetMatrixArray();
689 AMultB(ap,na,ncolsa,bp,nb,ncolsb,cp);
697template<
class Element>
703 if (
a.GetNcols() !=
b.GetNrows() ||
a.GetColLwb() !=
b.GetRowLwb()) {
704 Error(
"Mult",
"A rows and B columns incompatible");
708 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
709 Error(
"Mult",
"this->GetMatrixArray() == a.GetMatrixArray()");
713 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
714 Error(
"Mult",
"this->GetMatrixArray() == b.GetMatrixArray()");
720 const Element *ap =
a.GetMatrixArray();
721 const Element *bp =
b.GetMatrixArray();
722 Element *cp = this->GetMatrixArray();
723 if (
typeid(Element) ==
typeid(
Double_t))
724 cblas_dsymm (CblasRowMajor,CblasLeft,CblasUpper,fNrows,fNcols,1.0,
725 ap,
a.GetNcols(),bp,
b.GetNcols(),0.0,cp,fNcols);
726 else if (
typeid(Element) !=
typeid(
Float_t))
727 cblas_ssymm (CblasRowMajor,CblasLeft,CblasUpper,fNrows,fNcols,1.0,
728 ap,
a.GetNcols(),bp,
b.GetNcols(),0.0,cp,fNcols);
730 Error(
"Mult",
"type %s not implemented in BLAS library",
typeid(Element));
732 const Int_t na =
a.GetNoElements();
733 const Int_t nb =
b.GetNoElements();
734 const Int_t ncolsa =
a.GetNcols();
735 const Int_t ncolsb =
b.GetNcols();
736 const Element *
const ap =
a.GetMatrixArray();
737 const Element *
const bp =
b.GetMatrixArray();
738 Element * cp = this->GetMatrixArray();
740 AMultB(ap,na,ncolsa,bp,nb,ncolsb,cp);
749template<
class Element>
755 if (
a.GetNcols() !=
b.GetNrows() ||
a.GetColLwb() !=
b.GetRowLwb()) {
756 Error(
"Mult",
"A rows and B columns incompatible");
760 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
761 Error(
"Mult",
"this->GetMatrixArray() == a.GetMatrixArray()");
765 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
766 Error(
"Mult",
"this->GetMatrixArray() == b.GetMatrixArray()");
772 const Element *ap =
a.GetMatrixArray();
773 const Element *bp =
b.GetMatrixArray();
774 Element *cp = this->GetMatrixArray();
775 if (
typeid(Element) ==
typeid(
Double_t))
776 cblas_dsymm (CblasRowMajor,CblasRight,CblasUpper,fNrows,fNcols,1.0,
777 bp,
b.GetNcols(),ap,
a.GetNcols(),0.0,cp,fNcols);
778 else if (
typeid(Element) !=
typeid(
Float_t))
779 cblas_ssymm (CblasRowMajor,CblasRight,CblasUpper,fNrows,fNcols,1.0,
780 bp,
b.GetNcols(),ap,
a.GetNcols(),0.0,cp,fNcols);
782 Error(
"Mult",
"type %s not implemented in BLAS library",
typeid(Element));
784 const Int_t na =
a.GetNoElements();
785 const Int_t nb =
b.GetNoElements();
786 const Int_t ncolsa =
a.GetNcols();
787 const Int_t ncolsb =
b.GetNcols();
788 const Element *
const ap =
a.GetMatrixArray();
789 const Element *
const bp =
b.GetMatrixArray();
790 Element * cp = this->GetMatrixArray();
792 AMultB(ap,na,ncolsa,bp,nb,ncolsb,cp);
801template<
class Element>
807 if (
a.GetNcols() !=
b.GetNrows() ||
a.GetColLwb() !=
b.GetRowLwb()) {
808 Error(
"Mult",
"A rows and B columns incompatible");
812 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
813 Error(
"Mult",
"this->GetMatrixArray() == a.GetMatrixArray()");
817 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
818 Error(
"Mult",
"this->GetMatrixArray() == b.GetMatrixArray()");
824 const Element *ap =
a.GetMatrixArray();
825 const Element *bp =
b.GetMatrixArray();
826 Element *cp = this->GetMatrixArray();
827 if (
typeid(Element) ==
typeid(
Double_t))
828 cblas_dsymm (CblasRowMajor,CblasLeft,CblasUpper,fNrows,fNcols,1.0,
829 ap,
a.GetNcols(),bp,
b.GetNcols(),0.0,cp,fNcols);
830 else if (
typeid(Element) !=
typeid(
Float_t))
831 cblas_ssymm (CblasRowMajor,CblasLeft,CblasUpper,fNrows,fNcols,1.0,
832 ap,
a.GetNcols(),bp,
b.GetNcols(),0.0,cp,fNcols);
834 Error(
"Mult",
"type %s not implemented in BLAS library",
typeid(Element));
836 const Int_t na =
a.GetNoElements();
837 const Int_t nb =
b.GetNoElements();
838 const Int_t ncolsa =
a.GetNcols();
839 const Int_t ncolsb =
b.GetNcols();
840 const Element *
const ap =
a.GetMatrixArray();
841 const Element *
const bp =
b.GetMatrixArray();
842 Element * cp = this->GetMatrixArray();
844 AMultB(ap,na,ncolsa,bp,nb,ncolsb,cp);
852template<
class Element>
858 if (
a.GetNrows() !=
b.GetNrows() ||
a.GetRowLwb() !=
b.GetRowLwb()) {
859 Error(
"TMult",
"A rows and B columns incompatible");
863 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
864 Error(
"TMult",
"this->GetMatrixArray() == a.GetMatrixArray()");
868 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
869 Error(
"TMult",
"this->GetMatrixArray() == b.GetMatrixArray()");
875 const Element *ap =
a.GetMatrixArray();
876 const Element *bp =
b.GetMatrixArray();
877 Element *cp = this->GetMatrixArray();
878 if (
typeid(Element) ==
typeid(
Double_t))
879 cblas_dgemm (CblasRowMajor,CblasTrans,CblasNoTrans,this->fNrows,this->fNcols,
a.GetNrows(),
880 1.0,ap,
a.GetNcols(),bp,
b.GetNcols(),1.0,cp,this->fNcols);
881 else if (
typeid(Element) !=
typeid(
Float_t))
882 cblas_sgemm (CblasRowMajor,CblasTrans,CblasNoTrans,fNrows,fNcols,
a.GetNrows(),
883 1.0,ap,
a.GetNcols(),bp,
b.GetNcols(),1.0,cp,fNcols);
885 Error(
"TMult",
"type %s not implemented in BLAS library",
typeid(Element));
887 const Int_t nb =
b.GetNoElements();
888 const Int_t ncolsa =
a.GetNcols();
889 const Int_t ncolsb =
b.GetNcols();
890 const Element *
const ap =
a.GetMatrixArray();
891 const Element *
const bp =
b.GetMatrixArray();
892 Element * cp = this->GetMatrixArray();
894 AtMultB(ap,ncolsa,bp,nb,ncolsb,cp);
902template<
class Element>
908 if (
a.GetNrows() !=
b.GetNrows() ||
a.GetRowLwb() !=
b.GetRowLwb()) {
909 Error(
"TMult",
"A rows and B columns incompatible");
913 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
914 Error(
"TMult",
"this->GetMatrixArray() == a.GetMatrixArray()");
918 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
919 Error(
"TMult",
"this->GetMatrixArray() == b.GetMatrixArray()");
925 const Element *ap =
a.GetMatrixArray();
926 const Element *bp =
b.GetMatrixArray();
927 Element *cp = this->GetMatrixArray();
928 if (
typeid(Element) ==
typeid(
Double_t))
929 cblas_dgemm (CblasRowMajor,CblasTrans,CblasNoTrans,fNrows,fNcols,
a.GetNrows(),
930 1.0,ap,
a.GetNcols(),bp,
b.GetNcols(),1.0,cp,fNcols);
931 else if (
typeid(Element) !=
typeid(
Float_t))
932 cblas_sgemm (CblasRowMajor,CblasTrans,CblasNoTrans,fNrows,fNcols,
a.GetNrows(),
933 1.0,ap,
a.GetNcols(),bp,
b.GetNcols(),1.0,cp,fNcols);
935 Error(
"TMult",
"type %s not implemented in BLAS library",
typeid(Element));
937 const Int_t nb =
b.GetNoElements();
938 const Int_t ncolsa =
a.GetNcols();
939 const Int_t ncolsb =
b.GetNcols();
940 const Element *
const ap =
a.GetMatrixArray();
941 const Element *
const bp =
b.GetMatrixArray();
942 Element * cp = this->GetMatrixArray();
944 AtMultB(ap,ncolsa,bp,nb,ncolsb,cp);
951template<
class Element>
958 if (
a.GetNcols() !=
b.GetNcols() ||
a.GetColLwb() !=
b.GetColLwb()) {
959 Error(
"MultT",
"A rows and B columns incompatible");
963 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
964 Error(
"MultT",
"this->GetMatrixArray() == a.GetMatrixArray()");
968 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
969 Error(
"MultT",
"this->GetMatrixArray() == b.GetMatrixArray()");
975 const Element *ap =
a.GetMatrixArray();
976 const Element *bp =
b.GetMatrixArray();
977 Element *cp = this->GetMatrixArray();
978 if (
typeid(Element) ==
typeid(
Double_t))
979 cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasTrans,fNrows,fNcols,
a.GetNcols(),
980 1.0,ap,
a.GetNcols(),bp,
b.GetNcols(),1.0,cp,fNcols);
981 else if (
typeid(Element) !=
typeid(
Float_t))
982 cblas_sgemm (CblasRowMajor,CblasNoTrans,CblasTrans,fNrows,fNcols,
a.GetNcols(),
983 1.0,ap,
a.GetNcols(),bp,
b.GetNcols(),1.0,cp,fNcols);
985 Error(
"MultT",
"type %s not implemented in BLAS library",
typeid(Element));
987 const Int_t na =
a.GetNoElements();
988 const Int_t nb =
b.GetNoElements();
989 const Int_t ncolsa =
a.GetNcols();
990 const Int_t ncolsb =
b.GetNcols();
991 const Element *
const ap =
a.GetMatrixArray();
992 const Element *
const bp =
b.GetMatrixArray();
993 Element * cp = this->GetMatrixArray();
995 AMultBt(ap,na,ncolsa,bp,nb,ncolsb,cp);
1003template<
class Element>
1009 if (
a.GetNcols() !=
b.GetNcols() ||
a.GetColLwb() !=
b.GetColLwb()) {
1010 Error(
"MultT",
"A rows and B columns incompatible");
1014 if (this->GetMatrixArray() ==
a.GetMatrixArray()) {
1015 Error(
"MultT",
"this->GetMatrixArray() == a.GetMatrixArray()");
1019 if (this->GetMatrixArray() ==
b.GetMatrixArray()) {
1020 Error(
"MultT",
"this->GetMatrixArray() == b.GetMatrixArray()");
1026 const Element *ap =
a.GetMatrixArray();
1027 const Element *bp =
b.GetMatrixArray();
1028 Element *cp = this->GetMatrixArray();
1029 if (
typeid(Element) ==
typeid(
Double_t))
1030 cblas_dgemm (CblasRowMajor,CblasNoTrans,CblasTrans,this->fNrows,this->fNcols,
a.GetNcols(),
1031 1.0,ap,
a.GetNcols(),bp,
b.GetNcols(),1.0,cp,this->fNcols);
1032 else if (
typeid(Element) !=
typeid(
Float_t))
1033 cblas_sgemm (CblasRowMajor,CblasNoTrans,CblasTrans,fNrows,fNcols,
a.GetNcols(),
1034 1.0,ap,
a.GetNcols(),bp,
b.GetNcols(),1.0,cp,fNcols);
1036 Error(
"MultT",
"type %s not implemented in BLAS library",
typeid(Element));
1038 const Int_t na =
a.GetNoElements();
1039 const Int_t nb =
b.GetNoElements();
1040 const Int_t ncolsa =
a.GetNcols();
1041 const Int_t ncolsb =
b.GetNcols();
1042 const Element *
const ap =
a.GetMatrixArray();
1043 const Element *
const bp =
b.GetMatrixArray();
1044 Element * cp = this->GetMatrixArray();
1046 AMultBt(ap,na,ncolsa,bp,nb,ncolsb,cp);
1053template<
class Element>
1058 if (row_upb < row_lwb)
1060 Error(
"Use",
"row_upb=%d < row_lwb=%d",row_upb,row_lwb);
1063 if (col_upb < col_lwb)
1065 Error(
"Use",
"col_upb=%d < col_lwb=%d",col_upb,col_lwb);
1071 this->fNrows = row_upb-row_lwb+1;
1072 this->fNcols = col_upb-col_lwb+1;
1073 this->fRowLwb = row_lwb;
1074 this->fColLwb = col_lwb;
1075 this->fNelems = this->fNrows*this->fNcols;
1089template<
class Element>
1095 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
1096 Error(
"GetSub",
"row_lwb out of bounds");
1099 if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
1100 Error(
"GetSub",
"col_lwb out of bounds");
1103 if (row_upb < this->fRowLwb || row_upb > this->fRowLwb+this->fNrows-1) {
1104 Error(
"GetSub",
"row_upb out of bounds");
1107 if (col_upb < this->fColLwb || col_upb > this->fColLwb+this->fNcols-1) {
1108 Error(
"GetSub",
"col_upb out of bounds");
1111 if (row_upb < row_lwb || col_upb < col_lwb) {
1112 Error(
"GetSub",
"row_upb < row_lwb || col_upb < col_lwb");
1121 const Int_t row_lwb_sub = (shift) ? 0 : row_lwb;
1122 const Int_t row_upb_sub = (shift) ? row_upb-row_lwb : row_upb;
1123 const Int_t col_lwb_sub = (shift) ? 0 : col_lwb;
1124 const Int_t col_upb_sub = (shift) ? col_upb-col_lwb : col_upb;
1126 target.
ResizeTo(row_lwb_sub,row_upb_sub,col_lwb_sub,col_upb_sub);
1127 const Int_t nrows_sub = row_upb_sub-row_lwb_sub+1;
1128 const Int_t ncols_sub = col_upb_sub-col_lwb_sub+1;
1131 for (
Int_t irow = 0; irow < nrows_sub; irow++) {
1132 for (
Int_t icol = 0; icol < ncols_sub; icol++) {
1133 target(irow+row_lwb_sub,icol+col_lwb_sub) = (*this)(row_lwb+irow,col_lwb+icol);
1137 const Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNcols+(col_lwb-this->fColLwb);
1140 for (
Int_t irow = 0; irow < nrows_sub; irow++) {
1141 const Element *ap_sub = ap;
1142 for (
Int_t icol = 0; icol < ncols_sub; icol++) {
1156template<
class Element>
1163 if (row_lwb < this->fRowLwb || row_lwb > this->fRowLwb+this->fNrows-1) {
1164 Error(
"SetSub",
"row_lwb outof bounds");
1167 if (col_lwb < this->fColLwb || col_lwb > this->fColLwb+this->fNcols-1) {
1168 Error(
"SetSub",
"col_lwb outof bounds");
1171 if (row_lwb+source.
GetNrows() > this->fRowLwb+this->fNrows ||
1172 col_lwb+source.
GetNcols() > this->fColLwb+this->fNcols) {
1173 Error(
"SetSub",
"source matrix too large");
1184 for (
Int_t irow = 0; irow < nRows_source; irow++) {
1185 for (
Int_t icol = 0; icol < nCols_source; icol++) {
1186 (*this)(row_lwb+irow,col_lwb+icol) = source(rowlwb_s+irow,collwb_s+icol);
1191 Element *ap = this->GetMatrixArray()+(row_lwb-this->fRowLwb)*this->fNcols+(col_lwb-this->fColLwb);
1193 for (
Int_t irow = 0; irow < nRows_source; irow++) {
1194 Element *ap_sub = ap;
1195 for (
Int_t icol = 0; icol < nCols_source; icol++) {
1210template<
class Element>
1214 if (!this->fIsOwner) {
1215 Error(
"ResizeTo(Int_t,Int_t)",
"Not owner of data array,cannot resize");
1219 if (this->fNelems > 0) {
1220 if (this->fNrows == nrows && this->fNcols == ncols)
1222 else if (nrows == 0 || ncols == 0) {
1223 this->fNrows = nrows; this->fNcols = ncols;
1228 Element *elements_old = GetMatrixArray();
1229 const Int_t nelems_old = this->fNelems;
1230 const Int_t nrows_old = this->fNrows;
1231 const Int_t ncols_old = this->fNcols;
1233 Allocate(nrows,ncols);
1236 Element *elements_new = GetMatrixArray();
1239 if (this->fNelems > this->kSizeMax || nelems_old > this->kSizeMax)
1240 memset(elements_new,0,this->fNelems*
sizeof(Element));
1241 else if (this->fNelems > nelems_old)
1242 memset(elements_new+nelems_old,0,(this->fNelems-nelems_old)*
sizeof(Element));
1248 const Int_t nelems_new = this->fNelems;
1249 if (ncols_old < this->fNcols) {
1250 for (
Int_t i = nrows_copy-1; i >= 0; i--) {
1251 Memcpy_m(elements_new+i*this->fNcols,elements_old+i*ncols_old,ncols_copy,
1252 nelems_new,nelems_old);
1253 if (this->fNelems <= this->kSizeMax && nelems_old <= this->kSizeMax)
1254 memset(elements_new+i*this->fNcols+ncols_copy,0,(this->fNcols-ncols_copy)*
sizeof(Element));
1257 for (
Int_t i = 0; i < nrows_copy; i++)
1258 Memcpy_m(elements_new+i*this->fNcols,elements_old+i*ncols_old,ncols_copy,
1259 nelems_new,nelems_old);
1262 Delete_m(nelems_old,elements_old);
1264 Allocate(nrows,ncols,0,0,1);
1275template<
class Element>
1280 if (!this->fIsOwner) {
1281 Error(
"ResizeTo(Int_t,Int_t,Int_t,Int_t)",
"Not owner of data array,cannot resize");
1285 const Int_t new_nrows = row_upb-row_lwb+1;
1286 const Int_t new_ncols = col_upb-col_lwb+1;
1288 if (this->fNelems > 0) {
1290 if (this->fNrows == new_nrows && this->fNcols == new_ncols &&
1291 this->fRowLwb == row_lwb && this->fColLwb == col_lwb)
1293 else if (new_nrows == 0 || new_ncols == 0) {
1294 this->fNrows = new_nrows; this->fNcols = new_ncols;
1295 this->fRowLwb = row_lwb; this->fColLwb = col_lwb;
1300 Element *elements_old = GetMatrixArray();
1301 const Int_t nelems_old = this->fNelems;
1302 const Int_t nrows_old = this->fNrows;
1303 const Int_t ncols_old = this->fNcols;
1304 const Int_t rowLwb_old = this->fRowLwb;
1305 const Int_t colLwb_old = this->fColLwb;
1307 Allocate(new_nrows,new_ncols,row_lwb,col_lwb);
1310 Element *elements_new = GetMatrixArray();
1313 if (this->fNelems > this->kSizeMax || nelems_old > this->kSizeMax)
1314 memset(elements_new,0,this->fNelems*
sizeof(Element));
1315 else if (this->fNelems > nelems_old)
1316 memset(elements_new+nelems_old,0,(this->fNelems-nelems_old)*
sizeof(Element));
1321 const Int_t rowUpb_copy =
TMath::Min(this->fRowLwb+this->fNrows-1,rowLwb_old+nrows_old-1);
1322 const Int_t colUpb_copy =
TMath::Min(this->fColLwb+this->fNcols-1,colLwb_old+ncols_old-1);
1324 const Int_t nrows_copy = rowUpb_copy-rowLwb_copy+1;
1325 const Int_t ncols_copy = colUpb_copy-colLwb_copy+1;
1327 if (nrows_copy > 0 && ncols_copy > 0) {
1328 const Int_t colOldOff = colLwb_copy-colLwb_old;
1329 const Int_t colNewOff = colLwb_copy-this->fColLwb;
1330 if (ncols_old < this->fNcols) {
1331 for (
Int_t i = nrows_copy-1; i >= 0; i--) {
1332 const Int_t iRowOld = rowLwb_copy+i-rowLwb_old;
1333 const Int_t iRowNew = rowLwb_copy+i-this->fRowLwb;
1334 Memcpy_m(elements_new+iRowNew*this->fNcols+colNewOff,
1335 elements_old+iRowOld*ncols_old+colOldOff,ncols_copy,this->fNelems,nelems_old);
1336 if (this->fNelems <= this->kSizeMax && nelems_old <= this->kSizeMax)
1337 memset(elements_new+iRowNew*this->fNcols+colNewOff+ncols_copy,0,
1338 (this->fNcols-ncols_copy)*
sizeof(Element));
1341 for (
Int_t i = 0; i < nrows_copy; i++) {
1342 const Int_t iRowOld = rowLwb_copy+i-rowLwb_old;
1343 const Int_t iRowNew = rowLwb_copy+i-this->fRowLwb;
1344 Memcpy_m(elements_new+iRowNew*this->fNcols+colNewOff,
1345 elements_old+iRowOld*ncols_old+colOldOff,ncols_copy,this->fNelems,nelems_old);
1350 Delete_m(nelems_old,elements_old);
1352 Allocate(new_nrows,new_ncols,row_lwb,col_lwb,1);
1361template<
class Element>
1374template<
class Element>
1396template<
class Element>
1410template<
class Element>
1419 if (this->GetNrows() != this->GetNcols() || this->GetRowLwb() != this->GetColLwb()) {
1420 Error(
"Invert()",
"matrix should be square");
1422 Element *pM = this->GetMatrixArray();
1424 Error(
"InvertFast",
"matrix is singular");
1436 TMatrixTCramerInv::Inv2x2<Element>(*
this,det);
1441 TMatrixTCramerInv::Inv3x3<Element>(*
this,det);
1446 TMatrixTCramerInv::Inv4x4<Element>(*
this,det);
1451 TMatrixTCramerInv::Inv5x5<Element>(*
this,det);
1456 TMatrixTCramerInv::Inv6x6<Element>(*
this,det);
1469template<
class Element>
1476 Element *ap = this->GetMatrixArray();
1477 if (this->fNrows == this->fNcols && this->fRowLwb == this->fColLwb) {
1478 for (
Int_t i = 0; i < this->fNrows; i++) {
1479 const Int_t off_i = i*this->fNrows;
1480 for (
Int_t j = i+1; j < this->fNcols; j++) {
1481 const Int_t off_j = j*this->fNcols;
1482 const Element tmp = ap[off_i+j];
1483 ap[off_i+j] = ap[off_j+i];
1490 const Int_t nrows_old = this->fNrows;
1491 const Int_t ncols_old = this->fNcols;
1492 const Int_t rowlwb_old = this->fRowLwb;
1493 const Int_t collwb_old = this->fColLwb;
1495 this->fNrows = ncols_old; this->fNcols = nrows_old;
1496 this->fRowLwb = collwb_old; this->fColLwb = rowlwb_old;
1497 for (
Int_t irow = this->fRowLwb; irow < this->fRowLwb+this->fNrows; irow++) {
1498 for (
Int_t icol = this->fColLwb; icol < this->fColLwb+this->fNcols; icol++) {
1499 const Int_t off = (icol-collwb_old)*ncols_old;
1500 (*this)(irow,icol) = oldElems[off+irow-rowlwb_old];
1506 if (this->fNrows != source.
GetNcols() || this->fNcols != source.
GetNrows() ||
1509 Error(
"Transpose",
"matrix has wrong shape");
1514 const Element *scp = sp1;
1515 Element *tp = this->GetMatrixArray();
1516 const Element *
const tp_last = this->GetMatrixArray()+this->fNelems;
1520 while (tp < tp_last) {
1521 const Element *sp2 = scp++;
1524 while (sp2 < sp1+this->fNelems) {
1526 sp2 += this->fNrows;
1529 R__ASSERT(tp == tp_last && scp == sp1+this->fNrows);
1539template<
class Element>
1545 if (
v.GetNoElements() <
TMath::Max(this->fNrows,this->fNcols)) {
1546 Error(
"Rank1Update",
"vector too short");
1551 const Element *
const pv =
v.GetMatrixArray();
1552 Element *mp = this->GetMatrixArray();
1554 for (
Int_t i = 0; i < this->fNrows; i++) {
1555 const Element tmp = alpha*pv[i];
1556 for (
Int_t j = 0; j < this->fNcols; j++)
1567template<
class Element>
1574 if (
v1.GetNoElements() < this->fNrows) {
1575 Error(
"Rank1Update",
"vector v1 too short");
1579 if (
v2.GetNoElements() < this->fNcols) {
1580 Error(
"Rank1Update",
"vector v2 too short");
1585 const Element *
const pv1 =
v1.GetMatrixArray();
1586 const Element *
const pv2 =
v2.GetMatrixArray();
1587 Element *mp = this->GetMatrixArray();
1589 for (
Int_t i = 0; i < this->fNrows; i++) {
1590 const Element tmp = alpha*pv1[i];
1591 for (
Int_t j = 0; j < this->fNcols; j++)
1592 *mp++ += tmp*pv2[j];
1601template<
class Element>
1607 if (this->fNcols != this->fNrows || this->fColLwb != this->fRowLwb) {
1608 Error(
"Similarity(const TVectorT &)",
"matrix is not square");
1612 if (this->fNcols !=
v.GetNrows() || this->fColLwb !=
v.GetLwb()) {
1613 Error(
"Similarity(const TVectorT &)",
"vector and matrix incompatible");
1618 const Element *mp = this->GetMatrixArray();
1619 const Element *vp =
v.GetMatrixArray();
1622 const Element *
const vp_first = vp;
1623 const Element *
const vp_last = vp+
v.GetNrows();
1624 while (vp < vp_last) {
1626 for (
const Element *sp = vp_first; sp < vp_last; )
1627 sum2 += *mp++ * *sp++;
1628 sum1 += sum2 * *vp++;
1631 R__ASSERT(mp == this->GetMatrixArray()+this->GetNoElements());
1642template<
class Element>
1648 if (
v.GetNoElements() < this->fNrows) {
1649 Error(
"NormByColumn",
"vector shorter than matrix column");
1658 const Element *pv =
v.GetMatrixArray();
1659 Element *mp = this->GetMatrixArray();
1660 const Element *
const mp_last = mp+this->fNelems;
1663 for ( ; mp < mp_last; pv++) {
1664 for (
Int_t j = 0; j < this->fNcols; j++)
1669 Error(
"NormbyColumn",
"vector element %ld is zero",
Long_t(pv-
v.GetMatrixArray()));
1675 for ( ; mp < mp_last; pv++)
1676 for (
Int_t j = 0; j < this->fNcols; j++)
1689template<
class Element>
1695 if (
v.GetNoElements() < this->fNcols) {
1696 Error(
"NormByRow",
"vector shorter than matrix column");
1705 const Element *pv0 =
v.GetMatrixArray();
1706 const Element *pv = pv0;
1707 Element *mp = this->GetMatrixArray();
1708 const Element *
const mp_last = mp+this->fNelems;
1711 for ( ; mp < mp_last; pv = pv0 )
1712 for (
Int_t j = 0; j < this->fNcols; j++) {
1716 Error(
"NormbyRow",
"vector element %ld is zero",
Long_t(pv-pv0));
1721 for ( ; mp < mp_last; pv = pv0 )
1722 for (
Int_t j = 0; j < this->fNcols; j++)
1732template<
class Element>
1736 Error(
"operator=(const TMatrixT &)",
"matrices not compatible");
1742 memcpy(fElements,source.
GetMatrixArray(),this->fNelems*
sizeof(Element));
1743 this->fTol = source.
GetTol();
1751template<
class Element>
1755 Error(
"operator=(const TMatrixTSym &)",
"matrices not compatible");
1761 memcpy(fElements,source.
GetMatrixArray(),this->fNelems*
sizeof(Element));
1762 this->fTol = source.
GetTol();
1770template<
class Element>
1774 this->GetNrows() != source.
GetNrows()) || this->GetNcols() != source.
GetNcols() ||
1775 this->GetRowLwb() != source.
GetRowLwb() || this->GetColLwb() != source.
GetColLwb()) {
1776 Error(
"operator=(const TMatrixTSparse &",
"matrices not compatible");
1782 memset(fElements,0,this->fNelems*
sizeof(Element));
1785 Element * tp = this->GetMatrixArray();
1790 for (
Int_t irow = 0; irow < this->fNrows; irow++ ) {
1791 const Int_t off = irow*this->fNcols;
1792 const Int_t sIndex = pRowIndex[irow];
1793 const Int_t eIndex = pRowIndex[irow+1];
1794 for (
Int_t index = sIndex; index < eIndex; index++)
1795 tp[off+pColIndex[index]] = sp[index];
1797 this->fTol = source.
GetTol();
1805template<
class Element>
1810 if (lazy_constructor.
GetRowUpb() != this->GetRowUpb() ||
1811 lazy_constructor.
GetColUpb() != this->GetColUpb() ||
1812 lazy_constructor.
GetRowLwb() != this->GetRowLwb() ||
1813 lazy_constructor.
GetColLwb() != this->GetColLwb()) {
1814 Error(
"operator=(const TMatrixTLazy&)",
"matrix is incompatible with "
1815 "the assigned Lazy matrix");
1819 lazy_constructor.
FillIn(*
this);
1826template<
class Element>
1831 Element *ep = this->GetMatrixArray();
1832 const Element *
const ep_last = ep+this->fNelems;
1833 while (ep < ep_last)
1842template<
class Element>
1847 Element *ep = this->GetMatrixArray();
1848 const Element *
const ep_last = ep+this->fNelems;
1849 while (ep < ep_last)
1858template<
class Element>
1863 Element *ep = this->GetMatrixArray();
1864 const Element *
const ep_last = ep+this->fNelems;
1865 while (ep < ep_last)
1874template<
class Element>
1879 Element *ep = this->GetMatrixArray();
1880 const Element *
const ep_last = ep+this->fNelems;
1881 while (ep < ep_last)
1890template<
class Element>
1894 Error(
"operator+=(const TMatrixT &)",
"matrices not compatible");
1899 Element *tp = this->GetMatrixArray();
1900 const Element *
const tp_last = tp+this->fNelems;
1901 while (tp < tp_last)
1910template<
class Element>
1914 Error(
"operator+=(const TMatrixTSym &)",
"matrices not compatible");
1919 Element *tp = this->GetMatrixArray();
1920 const Element *
const tp_last = tp+this->fNelems;
1921 while (tp < tp_last)
1930template<
class Element>
1934 Error(
"operator=-(const TMatrixT &)",
"matrices not compatible");
1939 Element *tp = this->GetMatrixArray();
1940 const Element *
const tp_last = tp+this->fNelems;
1941 while (tp < tp_last)
1950template<
class Element>
1954 Error(
"operator=-(const TMatrixTSym &)",
"matrices not compatible");
1959 Element *tp = this->GetMatrixArray();
1960 const Element *
const tp_last = tp+this->fNelems;
1961 while (tp < tp_last)
1972template<
class Element>
1980 Error(
"operator*=(const TMatrixT &)",
"source matrix has wrong shape");
1997 Element work[kWorkMax];
1999 Element *trp = work;
2000 if (this->fNcols > kWorkMax) {
2001 isAllocated =
kTRUE;
2002 trp =
new Element[this->fNcols];
2005 Element *cp = this->GetMatrixArray();
2006 const Element *trp0 = cp;
2007 const Element *
const trp0_last = trp0+this->fNelems;
2008 while (trp0 < trp0_last) {
2009 memcpy(trp,trp0,this->fNcols*
sizeof(Element));
2010 for (
const Element *scp = sp; scp < sp+this->fNcols; ) {
2013 for (
Int_t j = 0; j < this->fNcols; j++) {
2014 cij += trp[j] * *scp;
2015 scp += this->fNcols;
2020 trp0 += this->fNcols;
2024 R__ASSERT(cp == trp0_last && trp0 == trp0_last);
2035template<
class Element>
2042 Error(
"operator*=(const TMatrixTSym &)",
"source matrix has wrong shape");
2059 Element work[kWorkMax];
2061 Element *trp = work;
2062 if (this->fNcols > kWorkMax) {
2063 isAllocated =
kTRUE;
2064 trp =
new Element[this->fNcols];
2067 Element *cp = this->GetMatrixArray();
2068 const Element *trp0 = cp;
2069 const Element *
const trp0_last = trp0+this->fNelems;
2070 while (trp0 < trp0_last) {
2071 memcpy(trp,trp0,this->fNcols*
sizeof(Element));
2072 for (
const Element *scp = sp; scp < sp+this->fNcols; ) {
2075 for (
Int_t j = 0; j < this->fNcols; j++) {
2076 cij += trp[j] * *scp;
2077 scp += this->fNcols;
2082 trp0 += this->fNcols;
2086 R__ASSERT(cp == trp0_last && trp0 == trp0_last);
2097template<
class Element>
2104 Error(
"operator*=(const TMatrixTDiag_const &)",
"wrong diagonal length");
2109 Element *mp = this->GetMatrixArray();
2110 const Element *
const mp_last = mp+this->fNelems;
2112 while (mp < mp_last) {
2113 const Element *dp = diag.
GetPtr();
2114 for (
Int_t j = 0; j < this->fNcols; j++) {
2127template<
class Element>
2134 Error(
"operator/=(const TMatrixTDiag_const &)",
"wrong diagonal length");
2139 Element *mp = this->GetMatrixArray();
2140 const Element *
const mp_last = mp+this->fNelems;
2142 while (mp < mp_last) {
2143 const Element *dp = diag.
GetPtr();
2144 for (
Int_t j = 0; j < this->fNcols; j++) {
2148 Error(
"operator/=",
"%d-diagonal element is zero",j);
2162template<
class Element>
2170 if (this->fNrows != mt->
GetNrows()) {
2171 Error(
"operator*=(const TMatrixTColumn_const &)",
"wrong column length");
2177 Element *mp = this->GetMatrixArray();
2178 const Element *
const mp_last = mp+this->fNelems;
2179 const Element *cp = col.
GetPtr();
2181 while (mp < mp_last) {
2183 for (
Int_t j = 0; j < this->fNcols; j++)
2195template<
class Element>
2203 if (this->fNrows != mt->
GetNrows()) {
2204 Error(
"operator/=(const TMatrixTColumn_const &)",
"wrong column matrix");
2210 Element *mp = this->GetMatrixArray();
2211 const Element *
const mp_last = mp+this->fNelems;
2212 const Element *cp = col.
GetPtr();
2214 while (mp < mp_last) {
2217 for (
Int_t j = 0; j < this->fNcols; j++)
2221 Error(
"operator/=",
"%d-row of matrix column is zero",icol);
2234template<
class Element>
2242 if (this->fNcols != mt->
GetNcols()) {
2243 Error(
"operator*=(const TMatrixTRow_const &)",
"wrong row length");
2249 Element *mp = this->GetMatrixArray();
2250 const Element *
const mp_last = mp+this->fNelems;
2252 while (mp < mp_last) {
2253 const Element *rp = row.
GetPtr();
2254 for (
Int_t j = 0; j < this->fNcols; j++) {
2268template<
class Element>
2275 if (this->fNcols != mt->
GetNcols()) {
2276 Error(
"operator/=(const TMatrixTRow_const &)",
"wrong row length");
2281 Element *mp = this->GetMatrixArray();
2282 const Element *
const mp_last = mp+this->fNelems;
2284 while (mp < mp_last) {
2285 const Element *rp = row.
GetPtr();
2286 for (
Int_t j = 0; j < this->fNcols; j++) {
2291 Error(
"operator/=",
"%d-col of matrix row is zero",j);
2307template<
class Element>
2310 if (!this->IsSymmetric())
2311 Warning(
"EigenVectors(TVectorT &)",
"Only real part of eigen-values will be returned");
2313 eigenValues.
ResizeTo(this->fNrows);
2321template<
class Element>
2332template<
class Element>
2343template<
class Element>
2352template<
class Element>
2363template<
class Element>
2372template<
class Element>
2383template<
class Element>
2394template<
class Element>
2397 return Element(-1.0)*(
operator-(source2,source1));
2403template<
class Element>
2414template<
class Element>
2417 return Element(-1.0)*
operator-(source,val);
2423template<
class Element>
2434template<
class Element>
2443template<
class Element>
2453template<
class Element>
2463template<
class Element>
2473template<
class Element>
2483template<
class Element>
2489 Error(
"operator&&(const TMatrixT&,const TMatrixT&)",
"matrices not compatible");
2499 while (tp < tp_last)
2500 *tp++ = (*sp1++ != 0.0 && *sp2++ != 0.0);
2508template<
class Element>
2514 Error(
"operator&&(const TMatrixT&,const TMatrixTSym&)",
"matrices not compatible");
2524 while (tp < tp_last)
2525 *tp++ = (*sp1++ != 0.0 && *sp2++ != 0.0);
2533template<
class Element>
2542template<
class Element>
2548 Error(
"operator||(const TMatrixT&,const TMatrixT&)",
"matrices not compatible");
2558 while (tp < tp_last)
2559 *tp++ = (*sp1++ != 0.0 || *sp2++ != 0.0);
2567template<
class Element>
2573 Error(
"operator||(const TMatrixT&,const TMatrixTSym&)",
"matrices not compatible");
2583 while (tp < tp_last)
2584 *tp++ = (*sp1++ != 0.0 || *sp2++ != 0.0);
2592template<
class Element>
2601template<
class Element>
2607 Error(
"operator|(const TMatrixT&,const TMatrixT&)",
"matrices not compatible");
2617 while (tp < tp_last) {
2618 *tp++ = (*sp1) > (*sp2); sp1++; sp2++;
2627template<
class Element>
2633 Error(
"operator>(const TMatrixT&,const TMatrixTSym&)",
"matrices not compatible");
2643 while (tp < tp_last) {
2644 *tp++ = (*sp1) > (*sp2); sp1++; sp2++;
2653template<
class Element>
2662template<
class Element>
2668 Error(
"operator>=(const TMatrixT&,const TMatrixT&)",
"matrices not compatible");
2678 while (tp < tp_last) {
2679 *tp++ = (*sp1) >= (*sp2); sp1++; sp2++;
2688template<
class Element>
2694 Error(
"operator>=(const TMatrixT&,const TMatrixTSym&)",
"matrices not compatible");
2704 while (tp < tp_last) {
2705 *tp++ = (*sp1) >= (*sp2); sp1++; sp2++;
2714template<
class Element>
2723template<
class Element>
2729 Error(
"operator<=(const TMatrixT&,const TMatrixT&)",
"matrices not compatible");
2739 while (tp < tp_last) {
2740 *tp++ = (*sp1) <= (*sp2); sp1++; sp2++;
2749template<
class Element>
2755 Error(
"operator<=(const TMatrixT&,const TMatrixTSym&)",
"matrices not compatible");
2765 while (tp < tp_last) {
2766 *tp++ = (*sp1) <= (*sp2); sp1++; sp2++;
2775template<
class Element>
2784template<
class Element>
2790 Error(
"operator<(const TMatrixT&,const TMatrixT&)",
"matrices not compatible");
2798 while (tp < tp_last) {
2799 *tp++ = (*sp1) < (*sp2); sp1++; sp2++;
2808template<
class Element>
2814 Error(
"operator<(const TMatrixT&,const TMatrixTSym&)",
"matrices not compatible");
2824 while (tp < tp_last) {
2825 *tp++ = (*sp1) < (*sp2); sp1++; sp2++;
2834template<
class Element>
2843template<
class Element>
2849 Error(
"operator!=(const TMatrixT&,const TMatrixT&)",
"matrices not compatible");
2859 while (tp != tp_last) {
2860 *tp++ = (*sp1) != (*sp2); sp1++; sp2++;
2869template<
class Element>
2875 Error(
"operator!=(const TMatrixT&,const TMatrixTSym&)",
"matrices not compatible");
2885 while (tp != tp_last) {
2886 *tp++ = (*sp1) != (*sp2); sp1++; sp2++;
2895template<
class Element>
2905template<class Element>
2906TMatrixT<Element> operator!=(const TMatrixT<Element> &source1,Element val)
2908 TMatrixT<Element> target; target.ResizeTo(source1);
2910 const Element *sp = source1.GetMatrixArray();
2911 Element *tp = target.GetMatrixArray();
2912 const Element * const tp_last = tp+target.GetNoElements();
2913 while (tp != tp_last) {
2914 *tp++ = (*sp != val); sp++;
2923template<class Element>
2924TMatrixT<Element> operator!=(Element val,const TMatrixT<Element> &source1)
2926 return operator!=(source1,val);
2933template<
class Element>
2937 ::Error(
"Add(TMatrixT &,Element,const TMatrixT &)",
"matrices not compatible");
2946 *tp++ = scalar * (*sp++);
2947 }
else if (scalar == 1.) {
2952 *tp++ += scalar * (*sp++);
2961template<
class Element>
2965 ::Error(
"Add(TMatrixT &,Element,const TMatrixTSym &)",
"matrices not compatible");
2973 *tp++ += scalar * (*sp++);
2981template<
class Element>
2985 ::Error(
"ElementMult(TMatrixT &,const TMatrixT &)",
"matrices not compatible");
3001template<
class Element>
3005 ::Error(
"ElementMult(TMatrixT &,const TMatrixTSym &)",
"matrices not compatible");
3021template<
class Element>
3025 ::Error(
"ElementDiv(TMatrixT &,const TMatrixT &)",
"matrices not compatible");
3032 while ( tp < ftp ) {
3038 Error(
"ElementDiv",
"source (%d,%d) is zero",irow,icol);
3049template<
class Element>
3053 ::Error(
"ElementDiv(TMatrixT &,const TMatrixTSym &)",
"matrices not compatible");
3060 while ( tp < ftp ) {
3066 Error(
"ElementDiv",
"source (%d,%d) is zero",irow,icol);
3077template<
class Element>
3079 const Element *
const bp,
Int_t nb,
Int_t ncolsb,Element *cp)
3081 const Element *arp0 = ap;
3082 while (arp0 < ap+na) {
3083 for (
const Element *bcp = bp; bcp < bp+ncolsb; ) {
3084 const Element *arp = arp0;
3086 while (bcp < bp+nb) {
3087 cij += *arp++ * *bcp;
3100template<
class Element>
3102 const Element *
const bp,
Int_t nb,
Int_t ncolsb,Element *cp)
3104 const Element *acp0 = ap;
3105 while (acp0 < ap+ncolsa) {
3106 for (
const Element *bcp = bp; bcp < bp+ncolsb; ) {
3107 const Element *acp = acp0;
3109 while (bcp < bp+nb) {
3124template<
class Element>
3126 const Element *
const bp,
Int_t nb,
Int_t ncolsb,Element *cp)
3128 const Element *arp0 = ap;
3129 while (arp0 < ap+na) {
3130 const Element *brp0 = bp;
3131 while (brp0 < bp+nb) {
3132 const Element *arp = arp0;
3133 const Element *brp = brp0;
3135 while (brp < brp0+ncolsb)
3136 cij += *arp++ * *brp++;
3147template<
class Element>
3156 }
else if (R__v == 2) {
3158 TObject::Streamer(R__b);
3160 R__b >> this->fNrows;
3161 R__b >> this->fNcols;
3162 R__b >> this->fNelems;
3163 R__b >> this->fRowLwb;
3164 R__b >> this->fColLwb;
3168 if (this->fNelems > 0) {
3169 fElements =
new Element[this->fNelems];
3176 TObject::Streamer(R__b);
3178 R__b >> this->fNrows;
3179 R__b >> this->fNcols;
3180 R__b >> this->fRowLwb;
3181 R__b >> this->fColLwb;
3182 this->fNelems = R__b.
ReadArray(fElements);
3186 if (R__v <= 2 && fElements) {
3187 for (
Int_t i = 0; i < this->fNrows; i++) {
3188 const Int_t off_i = i*this->fNcols;
3189 for (
Int_t j = i; j < this->fNcols; j++) {
3190 const Int_t off_j = j*this->fNrows;
3191 const Element tmp = fElements[off_i+j];
3192 fElements[off_i+j] = fElements[off_j+i];
3193 fElements[off_j+i] = tmp;
3197 if (this->fNelems > 0 && this->fNelems <= this->kSizeMax) {
3199 memcpy(fDataStack,fElements,this->fNelems*
sizeof(Element));
3200 delete [] fElements;
3202 fElements = fDataStack;
3203 }
else if (this->fNelems < 0)
3261template void TMatrixTAutoloadOps::AMultB <Float_t>(
const Float_t *
const ap,
Int_t na,
Int_t ncolsa,
3263template void TMatrixTAutoloadOps::AtMultB<Float_t>(
const Float_t *
const ap,
Int_t ncolsa,
3265template void TMatrixTAutoloadOps::AMultBt<Float_t>(
const Float_t *
const ap,
Int_t na,
Int_t ncolsa,
3318template void TMatrixTAutoloadOps::AMultB <Double_t>(
const Double_t *
const ap,
Int_t na,
Int_t ncolsa,
3320template void TMatrixTAutoloadOps::AtMultB<Double_t>(
const Double_t *
const ap,
Int_t ncolsa,
3322template void TMatrixTAutoloadOps::AMultBt<Double_t>(
const Double_t *
const ap,
Int_t na,
Int_t ncolsa,
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define templateClassImp(name)
Bool_t operator<=(const TDatime &d1, const TDatime &d2)
Bool_t operator>(const TDatime &d1, const TDatime &d2)
Bool_t operator!=(const TDatime &d1, const TDatime &d2)
Bool_t operator>=(const TDatime &d1, const TDatime &d2)
Bool_t operator<(const TDatime &d1, const TDatime &d2)
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
R__EXTERN Int_t gMatrixCheck
TString operator+(const TString &s1, const TString &s2)
Use the special concatenation constructor.
TTime operator*(const TTime &t1, const TTime &t2)
TTime operator-(const TTime &t1, const TTime &t2)
Buffer base class used for serializing objects.
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
virtual Int_t ReadArray(Bool_t *&b)=0
virtual Int_t CheckByteCount(UInt_t startpos, UInt_t bcnt, const TClass *clss)=0
virtual void ReadFastArray(Bool_t *b, Int_t n)=0
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=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.
const TMatrixD & GetEigenVectors() const
const TVectorD & GetEigenValuesRe() 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 > & SetMatrixArray(const Element *data, Option_t *option="")
Copy array data to matrix .
const TMatrixTBase< Element > * GetMatrix() const
const Element * GetPtr() const
const Element * GetPtr() const
const TMatrixTBase< Element > * GetMatrix() const
Templates of Lazy Matrix classes.
virtual void FillIn(TMatrixT< Element > &m) const =0
const Element * GetPtr() const
const TMatrixTBase< Element > * GetMatrix() const
virtual const Int_t * GetRowIndexArray() const
virtual const Element * GetMatrixArray() const
virtual const Int_t * GetColIndexArray() const
virtual const Element * GetMatrixArray() const
TMatrixT< Element > & Rank1Update(const TVectorT< Element > &v, Element alpha=1.0)
Perform a rank 1 operation on matrix A: A += alpha * v * v^T.
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])
void Delete_m(Int_t size, Element *&)
Delete data pointer m, if it was assigned on the heap.
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 > & operator*=(Element val)
Multiply every element of the matrix with val.
virtual const Element * GetMatrixArray() const
TMatrixT< Element > & operator=(const TMatrixT< Element > &source)
Assignment operator.
TMatrixT< Element > & NormByRow(const TVectorT< Element > &v, Option_t *option="D")
Multiply/divide matrix rows with a vector: option: "D" : b(i,j) = a(i,j)/v(j) i = 0,...
void Minus(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
General matrix summation. Create a matrix C such that C = A - B.
Element Similarity(const TVectorT< Element > &v) const
Calculate scalar v * (*this) * v^T.
TMatrixT< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant.
virtual Double_t Determinant() const
Return the matrix determinant.
void Plus(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
General matrix summation. Create a matrix C such that C = A + B.
TMatrixT< Element > & InvertFast(Double_t *det=0)
Invert the matrix and calculate its determinant, however upto (6x6) a fast Cramer inversion is used .
TMatrixT< Element > & operator-=(Element val)
Subtract val from every element of the matrix.
TMatrixT< Element > & Transpose(const TMatrixT< Element > &source)
Transpose matrix source.
void MultT(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
General matrix multiplication. Create a matrix C such that C = A * B^T.
virtual TMatrixTBase< Element > & SetSub(Int_t row_lwb, Int_t col_lwb, const TMatrixTBase< Element > &source)
Insert matrix source starting at [row_lwb][col_lwb], thereby overwriting the part [row_lwb....
virtual TMatrixTBase< Element > & GetSub(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, TMatrixTBase< Element > &target, Option_t *option="S") const
Get submatrix [row_lwb..row_upb] x [col_lwb..col_upb]; The indexing range of the returned matrix depe...
Element * New_m(Int_t size)
Return data pointer .
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.
TMatrixT< Element > & operator+=(Element val)
Add val to every element of the matrix.
TMatrixT< Element > & NormByColumn(const TVectorT< Element > &v, Option_t *option="D")
Multiply/divide matrix columns by a vector: option: "D" : b(i,j) = a(i,j)/v(i) i = 0,...
void TMult(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
Create a matrix C such that C = A' * B.
const TMatrixT< Element > EigenVectors(TVectorT< Element > &eigenValues) const
Return a matrix containing the eigen-vectors ordered by descending values of Re^2+Im^2 of the complex...
void Mult(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
General matrix multiplication. Create a matrix C such that C = A * B.
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 .
TMatrixT< Element > & operator/=(const TMatrixTDiag_const< Element > &diag)
Divide a matrix row by the diagonal of another matrix matrix(i,j) /= diag(j)
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)
TMatrixT< Element > operator!=(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
logical operation source1 != source2
void AMultBt(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^T.
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.
TMatrixT< Element > operator+(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
operation this = source1+source2
void AtMultB(const Element *const ap, Int_t ncolsa, const Element *const bp, Int_t nb, Int_t ncolsb, Element *cp)
Elementary routine to calculate matrix multiplication A^T*B.
TMatrixT< Element > operator>=(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
logical operation source1 >= source2
TMatrixT< Element > & ElementMult(TMatrixT< Element > &target, const TMatrixT< Element > &source)
Multiply target by the source, element-by-element.
TMatrixT< Element > operator||(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
Logical OR.
TMatrixT< Element > operator<(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
logical operation source1 < source2
TMatrixT< Element > operator>(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
logical operation source1 > source2
TMatrixT< Element > & Add(TMatrixT< Element > &target, Element scalar, const TMatrixT< Element > &source)
Modify addition: target += scalar * source.
TMatrixT< Element > & ElementDiv(TMatrixT< Element > &target, const TMatrixT< Element > &source)
Divide target by the source, element-by-element.
Element1 Mult(const TVectorT< Element1 > &v1, const TMatrixT< Element2 > &m, const TVectorT< Element3 > &v2)
Perform v1 * M * v2, a scalar result.
TMatrixT< Element > operator-(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
operation this = source1-source2
TMatrixT< Element > operator<=(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
logical operation source1 <= source2
TMatrixT< Element > operator&&(const TMatrixT< Element > &source1, const TMatrixT< Element > &source2)
Logical AND.
TMatrixT< Element > operator*(Element val, const TMatrixT< Element > &source)
operation this = val*source
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 .