File indexing completed on 2026-09-16 09:27:14
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 #ifndef ROOT_TMatrixT
0013 #define ROOT_TMatrixT
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023 #include "TMatrixTBase.h"
0024 #include "TMatrixTUtils.h"
0025
0026 #ifdef CBLAS
0027 #include <vecLib/vBLAS.h>
0028
0029 #endif
0030
0031 #include "Rtypes.h"
0032 #include "TError.h"
0033
0034 #include <cassert>
0035
0036 template<class Element> class TMatrixTSym;
0037 template<class Element> class TMatrixTSparse;
0038 template<class Element> class TMatrixTLazy;
0039
0040 template<class Element> class TMatrixT : public TMatrixTBase<Element> {
0041
0042 protected:
0043
0044 Element fDataStack[TMatrixTBase<Element>::kSizeMax];
0045 Element *fElements;
0046
0047 Element *New_m (Int_t size);
0048 void Delete_m(Int_t size,Element*&);
0049 Int_t Memcpy_m(Element *newp,const Element *oldp,Int_t copySize,
0050 Int_t newSize,Int_t oldSize);
0051 void Allocate(Int_t nrows,Int_t ncols,Int_t row_lwb = 0,Int_t col_lwb = 0,Int_t init = 0,
0052 Int_t = -1);
0053
0054
0055 public:
0056
0057
0058 enum {kWorkMax = 100};
0059 enum EMatrixCreatorsOp1 { kZero,kUnit,kTransposed,kInverted,kAtA };
0060 enum EMatrixCreatorsOp2 { kMult,kTransposeMult,kInvMult,kMultTranspose,
0061
0062 #if defined(__clang__) && __clang_major__ < 20
0063 #pragma clang diagnostic push
0064 #pragma clang diagnostic ignored "-Wshadow"
0065 #endif
0066 kPlus,
0067 #if defined(__clang__) && __clang_major__ < 20
0068 #pragma clang diagnostic pop
0069 #endif
0070 kMinus };
0071
0072 TMatrixT(): fDataStack(), fElements(nullptr) { }
0073 TMatrixT(Int_t nrows,Int_t ncols);
0074 TMatrixT(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
0075 TMatrixT(Int_t nrows,Int_t ncols,const Element *data,Option_t *option="");
0076 TMatrixT(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,const Element *data,Option_t *option="");
0077 TMatrixT(const TMatrixT <Element> &another);
0078 TMatrixT(const TMatrixTSym <Element> &another);
0079 TMatrixT(const TMatrixTSparse<Element> &another);
0080 template <class Element2> TMatrixT(const TMatrixT<Element2> &another): fElements(nullptr)
0081 {
0082 R__ASSERT(another.IsValid());
0083 Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb());
0084 *this = another;
0085 }
0086
0087 TMatrixT(EMatrixCreatorsOp1 op,const TMatrixT<Element> &prototype);
0088 TMatrixT(const TMatrixT <Element> &a,EMatrixCreatorsOp2 op,const TMatrixT <Element> &b);
0089 TMatrixT(const TMatrixT <Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSym<Element> &b);
0090 TMatrixT(const TMatrixTSym <Element> &a,EMatrixCreatorsOp2 op,const TMatrixT <Element> &b);
0091 TMatrixT(const TMatrixTSym <Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSym<Element> &b);
0092 TMatrixT(const TMatrixTLazy<Element> &lazy_constructor);
0093
0094 ~TMatrixT() override { TMatrixT::Clear(); }
0095
0096
0097
0098 void Plus (const TMatrixT <Element> &a,const TMatrixT <Element> &b);
0099 void Plus (const TMatrixT <Element> &a,const TMatrixTSym<Element> &b);
0100 void Plus (const TMatrixTSym<Element> &a,const TMatrixT <Element> &b) { Plus(b,a); }
0101
0102 void Minus(const TMatrixT <Element> &a,const TMatrixT <Element> &b);
0103 void Minus(const TMatrixT <Element> &a,const TMatrixTSym<Element> &b);
0104 void Minus(const TMatrixTSym<Element> &a,const TMatrixT <Element> &b) { Minus(b,a); }
0105
0106 void Mult (const TMatrixT <Element> &a,const TMatrixT <Element> &b);
0107 void Mult (const TMatrixT <Element> &a,const TMatrixTSym<Element> &b);
0108 void Mult (const TMatrixTSym<Element> &a,const TMatrixT <Element> &b);
0109 void Mult (const TMatrixTSym<Element> &a,const TMatrixTSym<Element> &b);
0110
0111 void TMult(const TMatrixT <Element> &a,const TMatrixT <Element> &b);
0112 void TMult(const TMatrixT <Element> &a,const TMatrixTSym<Element> &b);
0113 void TMult(const TMatrixTSym<Element> &a,const TMatrixT <Element> &b) { Mult(a,b); }
0114 void TMult(const TMatrixTSym<Element> &a,const TMatrixTSym<Element> &b) { Mult(a,b); }
0115
0116 void MultT(const TMatrixT <Element> &a,const TMatrixT <Element> &b);
0117 void MultT(const TMatrixT <Element> &a,const TMatrixTSym<Element> &b) { Mult(a,b); }
0118 void MultT(const TMatrixTSym<Element> &a,const TMatrixT <Element> &b);
0119 void MultT(const TMatrixTSym<Element> &a,const TMatrixTSym<Element> &b) { Mult(a,b); }
0120
0121 inline void SetElement(Int_t rown, Int_t coln, Element val);
0122
0123 const Element *GetMatrixArray () const override;
0124 Element *GetMatrixArray () override;
0125 const Int_t *GetRowIndexArray() const override { return nullptr; }
0126 Int_t *GetRowIndexArray() override { return nullptr; }
0127 const Int_t *GetColIndexArray() const override { return nullptr; }
0128 Int_t *GetColIndexArray() override { return nullptr; }
0129
0130 TMatrixTBase<Element> &SetRowIndexArray(Int_t * ) override { MayNotUse("SetRowIndexArray(Int_t *)"); return *this; }
0131 TMatrixTBase<Element> &SetColIndexArray(Int_t * ) override { MayNotUse("SetColIndexArray(Int_t *)"); return *this; }
0132
0133 void Clear(Option_t * ="") override
0134 {
0135 if (this->fIsOwner)
0136 Delete_m(this->fNelems, fElements);
0137 else
0138 fElements = nullptr;
0139 this->fNelems = 0;
0140 }
0141
0142 TMatrixT <Element> &Use (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Element *data);
0143 const TMatrixT <Element> &Use (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,const Element *data) const
0144 { return (const TMatrixT<Element>&)
0145 ((const_cast<TMatrixT<Element> *>(this))->Use(row_lwb,row_upb,col_lwb,col_upb, const_cast<Element *>(data))); }
0146 TMatrixT <Element> &Use (Int_t nrows,Int_t ncols,Element *data);
0147 const TMatrixT <Element> &Use (Int_t nrows,Int_t ncols,const Element *data) const;
0148 TMatrixT <Element> &Use (TMatrixT<Element> &a);
0149 const TMatrixT <Element> &Use (const TMatrixT<Element> &a) const;
0150
0151 TMatrixTBase<Element> &GetSub (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
0152 TMatrixTBase<Element> &target,Option_t *option="S") const override;
0153 TMatrixT <Element> GetSub (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Option_t *option="S") const;
0154 TMatrixTBase<Element> &SetSub (Int_t row_lwb,Int_t col_lwb,const TMatrixTBase<Element> &source) override;
0155
0156 TMatrixTBase<Element> &ResizeTo(Int_t nrows,Int_t ncols,Int_t =-1) override;
0157 TMatrixTBase<Element> &ResizeTo(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t =-1) override;
0158 inline TMatrixTBase<Element> &ResizeTo(const TMatrixT<Element> &m) {
0159 return ResizeTo(m.GetRowLwb(),m.GetRowUpb(),m.GetColLwb(),m.GetColUpb());
0160 }
0161
0162 Double_t Determinant () const override;
0163 void Determinant (Double_t &d1,Double_t &d2) const override;
0164
0165 TMatrixT<Element> &Invert (Double_t *det = nullptr);
0166 TMatrixT<Element> &InvertFast (Double_t *det = nullptr);
0167 TMatrixT<Element> &Transpose (const TMatrixT<Element> &source);
0168 inline TMatrixT<Element> &T () { return this->Transpose(*this); }
0169 TMatrixT<Element> &Rank1Update (const TVectorT<Element> &v,Element alpha=1.0);
0170 TMatrixT<Element> &Rank1Update (const TVectorT<Element> &v1,const TVectorT<Element> &v2,Element alpha=1.0);
0171 Element Similarity (const TVectorT<Element> &v) const;
0172
0173 TMatrixT<Element> &NormByColumn(const TVectorT<Element> &v,Option_t *option="D");
0174 TMatrixT<Element> &NormByRow (const TVectorT<Element> &v,Option_t *option="D");
0175
0176
0177 inline Element operator()(Int_t rown,Int_t coln) const override;
0178 inline Element &operator()(Int_t rown,Int_t coln) override;
0179
0180
0181 inline const TMatrixTRow_const<Element> operator[](Int_t rown) const { return TMatrixTRow_const<Element>(*this,rown); }
0182 inline TMatrixTRow <Element> operator[](Int_t rown) { return TMatrixTRow <Element>(*this,rown); }
0183
0184 TMatrixT<Element> &operator= (const TMatrixT <Element> &source);
0185 TMatrixT<Element> &operator= (const TMatrixTSym <Element> &source);
0186 TMatrixT<Element> &operator= (const TMatrixTSparse<Element> &source);
0187 TMatrixT<Element> &operator= (const TMatrixTLazy <Element> &source);
0188 template <class Element2> TMatrixT<Element> &operator= (const TMatrixT<Element2> &source)
0189 {
0190 if (!AreCompatible(*this,source)) {
0191 Error("operator=(const TMatrixT2 &)","matrices not compatible");
0192 return *this;
0193 }
0194
0195 TObject::operator=(source);
0196 const Element2 * const ps = source.GetMatrixArray();
0197 Element * const pt = this->GetMatrixArray();
0198 for (Int_t i = 0; i < this->fNelems; i++)
0199 pt[i] = ps[i];
0200 this->fTol = source.GetTol();
0201 return *this;
0202 }
0203
0204 TMatrixT<Element> &operator= (Element val);
0205 TMatrixT<Element> &operator-=(Element val);
0206 TMatrixT<Element> &operator+=(Element val);
0207 TMatrixT<Element> &operator*=(Element val);
0208
0209 TMatrixT<Element> &operator+=(const TMatrixT <Element> &source);
0210 TMatrixT<Element> &operator+=(const TMatrixTSym<Element> &source);
0211 TMatrixT<Element> &operator-=(const TMatrixT <Element> &source);
0212 TMatrixT<Element> &operator-=(const TMatrixTSym<Element> &source);
0213
0214 TMatrixT<Element> &operator*=(const TMatrixT <Element> &source);
0215 TMatrixT<Element> &operator*=(const TMatrixTSym <Element> &source);
0216 TMatrixT<Element> &operator*=(const TMatrixTDiag_const <Element> &diag);
0217 TMatrixT<Element> &operator/=(const TMatrixTDiag_const <Element> &diag);
0218 TMatrixT<Element> &operator*=(const TMatrixTRow_const <Element> &row);
0219 TMatrixT<Element> &operator/=(const TMatrixTRow_const <Element> &row);
0220 TMatrixT<Element> &operator*=(const TMatrixTColumn_const<Element> &col);
0221 TMatrixT<Element> &operator/=(const TMatrixTColumn_const<Element> &col);
0222
0223 const TMatrixT<Element> EigenVectors(TVectorT<Element> &eigenValues) const;
0224
0225 ClassDefOverride(TMatrixT,4)
0226 };
0227
0228
0229
0230
0231
0232
0233
0234
0235 template <> TClass *TMatrixT<double>::Class();
0236
0237
0238 template <class Element> inline const Element *TMatrixT<Element>::GetMatrixArray() const { return fElements; }
0239 template <class Element> inline Element *TMatrixT<Element>::GetMatrixArray() { return fElements; }
0240
0241 template <class Element> inline TMatrixT<Element> &TMatrixT<Element>::Use (Int_t nrows,Int_t ncols,Element *data)
0242 { return Use(0,nrows-1,0,ncols-1,data); }
0243 template <class Element> inline const TMatrixT<Element> &TMatrixT<Element>::Use (Int_t nrows,Int_t ncols,const Element *data) const
0244 { return Use(0,nrows-1,0,ncols-1,data); }
0245 template <class Element> inline TMatrixT<Element> &TMatrixT<Element>::Use (TMatrixT &a)
0246 {
0247 R__ASSERT(a.IsValid());
0248 return Use(a.GetRowLwb(),a.GetRowUpb(),
0249 a.GetColLwb(),a.GetColUpb(),a.GetMatrixArray());
0250 }
0251 template <class Element> inline const TMatrixT<Element> &TMatrixT<Element>::Use (const TMatrixT &a) const
0252 {
0253 R__ASSERT(a.IsValid());
0254 return Use(a.GetRowLwb(),a.GetRowUpb(),
0255 a.GetColLwb(),a.GetColUpb(),a.GetMatrixArray());
0256 }
0257
0258 template <class Element> inline TMatrixT<Element> TMatrixT<Element>::GetSub (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
0259 Option_t *option) const
0260 {
0261 TMatrixT tmp;
0262 this->GetSub(row_lwb,row_upb,col_lwb,col_upb,tmp,option);
0263 return tmp;
0264 }
0265
0266 template <class Element> inline Element TMatrixT<Element>::operator()(Int_t rown,Int_t coln) const
0267 {
0268 R__ASSERT(this->IsValid());
0269 const Int_t arown = rown-this->fRowLwb;
0270 const Int_t acoln = coln-this->fColLwb;
0271 if (arown >= this->fNrows || arown < 0) {
0272 Error("operator()","Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->fNrows);
0273 return TMatrixTBase<Element>::NaNValue();
0274 }
0275 if (acoln >= this->fNcols || acoln < 0) {
0276 Error("operator()","Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->fNcols);
0277 return TMatrixTBase<Element>::NaNValue();
0278
0279 }
0280 return (fElements[arown*this->fNcols+acoln]);
0281 }
0282
0283 template <class Element> inline Element &TMatrixT<Element>::operator()(Int_t rown,Int_t coln)
0284 {
0285 R__ASSERT(this->IsValid());
0286 const Int_t arown = rown-this->fRowLwb;
0287 const Int_t acoln = coln-this->fColLwb;
0288 if (arown >= this->fNrows || arown < 0) {
0289 Error("operator()","Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->fNrows);
0290 return TMatrixTBase<Element>::NaNValue();
0291 }
0292 if (acoln >= this->fNcols || acoln < 0) {
0293 Error("operator()","Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->fNcols);
0294 return TMatrixTBase<Element>::NaNValue();
0295 }
0296 return (fElements[arown*this->fNcols+acoln]);
0297 }
0298
0299
0300
0301
0302
0303 template <class Element>
0304 inline void TMatrixT<Element>::SetElement(Int_t rown, Int_t coln, Element val)
0305 {
0306 assert(this->IsValid());
0307 rown = rown - this->fRowLwb;
0308 coln = coln - this->fColLwb;
0309 assert((rown < this->fNrows && rown >= 0) && "SetElement() error: row index outside matrix range");
0310 assert((coln < this->fNcols && coln >= 0) && "SetElement() error: column index outside matrix range");
0311 fElements[rown * this->fNcols + coln] = val;
0312 }
0313
0314 inline namespace TMatrixTAutoloadOps {
0315
0316 template <class Element> TMatrixT<Element> operator+ (const TMatrixT <Element> &source1,const TMatrixT <Element> &source2);
0317 template <class Element> TMatrixT<Element> operator+ (const TMatrixT <Element> &source1,const TMatrixTSym<Element> &source2);
0318 template <class Element> TMatrixT<Element> operator+ (const TMatrixTSym<Element> &source1,const TMatrixT <Element> &source2);
0319 template <class Element> TMatrixT<Element> operator+ (const TMatrixT <Element> &source , Element val );
0320 template <class Element> TMatrixT<Element> operator+ ( Element val ,const TMatrixT <Element> &source );
0321 template <class Element> TMatrixT<Element> operator- (const TMatrixT <Element> &source1,const TMatrixT <Element> &source2);
0322 template <class Element> TMatrixT<Element> operator- (const TMatrixT <Element> &source1,const TMatrixTSym<Element> &source2);
0323 template <class Element> TMatrixT<Element> operator- (const TMatrixTSym<Element> &source1,const TMatrixT <Element> &source2);
0324 template <class Element> TMatrixT<Element> operator- (const TMatrixT <Element> &source , Element val );
0325 template <class Element> TMatrixT<Element> operator- ( Element val ,const TMatrixT <Element> &source );
0326 template <class Element> TMatrixT<Element> operator* ( Element val ,const TMatrixT <Element> &source );
0327 template <class Element> TMatrixT<Element> operator* (const TMatrixT <Element> &source , Element val );
0328 template <class Element> TMatrixT<Element> operator* (const TMatrixT <Element> &source1,const TMatrixT <Element> &source2);
0329 template <class Element> TMatrixT<Element> operator* (const TMatrixT <Element> &source1,const TMatrixTSym<Element> &source2);
0330 template <class Element> TMatrixT<Element> operator* (const TMatrixTSym<Element> &source1,const TMatrixT <Element> &source2);
0331 template <class Element> TMatrixT<Element> operator* (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element> &source2);
0332
0333 #if (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) >= 40600
0334 #pragma GCC diagnostic push
0335 #pragma GCC diagnostic ignored "-Weffc++"
0336 #endif
0337 template <class Element> TMatrixT<Element> operator&& (const TMatrixT <Element> &source1,const TMatrixT <Element> &source2);
0338 template <class Element> TMatrixT<Element> operator&& (const TMatrixT <Element> &source1,const TMatrixTSym<Element> &source2);
0339 template <class Element> TMatrixT<Element> operator&& (const TMatrixTSym<Element> &source1,const TMatrixT <Element> &source2);
0340 template <class Element> TMatrixT<Element> operator|| (const TMatrixT <Element> &source1,const TMatrixT <Element> &source2);
0341 template <class Element> TMatrixT<Element> operator|| (const TMatrixT <Element> &source1,const TMatrixTSym<Element> &source2);
0342 template <class Element> TMatrixT<Element> operator|| (const TMatrixTSym<Element> &source1,const TMatrixT <Element> &source2);
0343 #if (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) >= 40600
0344 #pragma GCC diagnostic pop
0345 #endif
0346 template <class Element> TMatrixT<Element> operator> (const TMatrixT <Element> &source1,const TMatrixT <Element> &source2);
0347 template <class Element> TMatrixT<Element> operator> (const TMatrixT <Element> &source1,const TMatrixTSym<Element> &source2);
0348 template <class Element> TMatrixT<Element> operator> (const TMatrixTSym<Element> &source1,const TMatrixT <Element> &source2);
0349 template <class Element> TMatrixT<Element> operator>= (const TMatrixT <Element> &source1,const TMatrixT <Element> &source2);
0350 template <class Element> TMatrixT<Element> operator>= (const TMatrixT <Element> &source1,const TMatrixTSym<Element> &source2);
0351 template <class Element> TMatrixT<Element> operator>= (const TMatrixTSym<Element> &source1,const TMatrixT <Element> &source2);
0352 template <class Element> TMatrixT<Element> operator<= (const TMatrixT <Element> &source1,const TMatrixT <Element> &source2);
0353 template <class Element> TMatrixT<Element> operator<= (const TMatrixT <Element> &source1,const TMatrixTSym<Element> &source2);
0354 template <class Element> TMatrixT<Element> operator<= (const TMatrixTSym<Element> &source1,const TMatrixT <Element> &source2);
0355 template <class Element> TMatrixT<Element> operator< (const TMatrixT <Element> &source1,const TMatrixT <Element> &source2);
0356 template <class Element> TMatrixT<Element> operator< (const TMatrixT <Element> &source1,const TMatrixTSym<Element> &source2);
0357 template <class Element> TMatrixT<Element> operator< (const TMatrixTSym<Element> &source1,const TMatrixT <Element> &source2);
0358 template <class Element> TMatrixT<Element> operator!= (const TMatrixT <Element> &source1,const TMatrixT <Element> &source2);
0359 template <class Element> TMatrixT<Element> operator!= (const TMatrixT <Element> &source1,const TMatrixTSym<Element> &source2);
0360 template <class Element> TMatrixT<Element> operator!= (const TMatrixTSym<Element> &source1,const TMatrixT <Element> &source2);
0361
0362 template <class Element> TMatrixT<Element> &Add (TMatrixT<Element> &target, Element scalar,const TMatrixT <Element> &source);
0363 template <class Element> TMatrixT<Element> &Add (TMatrixT<Element> &target, Element scalar,const TMatrixTSym<Element> &source);
0364 template <class Element> TMatrixT<Element> &ElementMult(TMatrixT<Element> &target,const TMatrixT <Element> &source);
0365 template <class Element> TMatrixT<Element> &ElementMult(TMatrixT<Element> &target,const TMatrixTSym<Element> &source);
0366 template <class Element> TMatrixT<Element> &ElementDiv (TMatrixT<Element> &target,const TMatrixT <Element> &source);
0367 template <class Element> TMatrixT<Element> &ElementDiv (TMatrixT<Element> &target,const TMatrixTSym<Element> &source);
0368
0369 template <class Element> void AMultB (const Element * const ap,Int_t na,Int_t ncolsa,
0370 const Element * const bp,Int_t nb,Int_t ncolsb,Element *cp);
0371 template <class Element> void AtMultB(const Element * const ap,Int_t ncolsa,
0372 const Element * const bp,Int_t nb,Int_t ncolsb,Element *cp);
0373 template <class Element> void AMultBt(const Element * const ap,Int_t na,Int_t ncolsa,
0374 const Element * const bp,Int_t nb,Int_t ncolsb,Element *cp);
0375 }
0376 #endif