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