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