TMatrixDBase : accessing fElements

In function VerifyMatrixIdentity, we need to access i,j value of m1. Through the below virtual call m1(i,j) we get access to the array fElement. Is there a direct way to get access to the array without making these virtual calls.

GetMatrixArray() can be used to access the array but it is still pure virtual in TMatrixDBase.

const Double_t dev = TMath::Abs(m1(i,j)-m2(i,j));

regards,
Hrishi

Hi Hrishi,

The signature of the VerifyMatrixIdentity is

Bool_t VerifyMatrixIdentity(const TMatrixDBase &m1,const TMatrixDBase &m2,…

It can be used to compare TMatrixD, TMatrixDSym and TMatrixDSparse
matrixes . Since they have different storage properties, we have
to access the matrix elements through m(i,j) or m[i][j] .

You seem to want to able to have access to the bare pointers
for speed ?? In that case, you will have to define your own
routine like for instance :

Bool_t VerifyMatrixIdentity(const TMatrixD &m1,const TMatrixD &m2,…

and replace

for (Int_t i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++) {
for (Int_t j = m1.GetColLwb(); j <= m1.GetColUpb(); j++) {
const Double_t dev = TMath::Abs(m1(i,j)-m2(i,j));
if (dev > maxDevObs) {
imax = i;
jmax = j;
maxDevObs = dev;
}
}
}

by

const Double_t * const p1 = m1.GetMatrixArray();
const Double_t * const p2 = m2.GetMatrixArray();
for (Int_t i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++) {
const Int_t off = i*m1.GetNrows();
for (Int_t j = m1.GetColLwb(); j <= m1.GetColUpb(); j++) {
const Double_t dev = TMath::Abs(p1[off+j]-p2[off+j]);
if (dev > maxDevObs) {
imax = i;
jmax = j;
maxDevObs = dev;
}
}
}

Eddy

Thanks for the help.

What I see is,

The implementation of GetMatrixArray is in TMatrixD, TMatrixDSym … It is still a pure virtual in TMatrixDBase.

What If I modify class TMatrixDBase so that I get fElements and GetMatrixArray implementation there itself.

example:

inline const Double_t *TMatrixDBase::GetMatrixArray() const { return fElements; }
inline Double_t *TMatrixDBase::GetMatrixArray() { return fElements; }

Thanks,
Hrishi

Hi,

Sure you could do that . However, the interpretation of the returned
data pointer is completely different for a TMatrixDSparse and
TMatrixD,TMatrixDSym . In the future, we could introduce a
difference between TMatrixD and TMatrixDSym .

Eddy

Thanks for the clarification.

I modified, TMatrixDBase.cxx

In verifymatrixidentity() added the following:

const Double_t * mp1 = m1.GetMatrixArray();
const Double_t * mp2 = m2.GetMatrixArray();
Int_t m1fRowLwb = m1.GetRowLwb();
Int_t m2fRowLwb = m2.GetRowLwb();
Int_t m1fColLwb = m1.GetColLwb();
Int_t m2fColLwb = m2.GetColLwb();
Int_t m1fNcols = m1.GetNcols ();
Int_t m2fNcols = m2.GetNcols ();

for (Int_t i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++) {
for (Int_t j = m1.GetColLwb(); j <= m1.GetColUpb(); j++) {
const Double_t dev = TMath::Abs(mp1[(i -m1fRowLwb)* m1fNcols+(j -m1fColLwb)]-mp2[(i -m2fRowLwb)* m2fNcols+(j -m2fColLwb)]);

and added the following to TMatrixDBase.h:

inline const Double_t *TMatrixDBase::GetMatrixArray() const { return fElements; }
inline Double_t *TMatrixDBase::GetMatrixArray() { return fElements; }

Compilation goes through fine, but get run time error in test, "StressLinear"
In functions mstress_matrix_io() and void astress_decomp_io(Int_t msize).

Donot know why we have this error in this function.

thanks,
Hrishi

Hi,

What you did does not make sense:

  1. You created now a routine that is in the base class and will
    not work for a Sparse matrix ??
  2. You changed a whole bunch of files and got an error .
    Do you now expect that others make the same changes and
    reproduce it ??

If you would have followed my suggestion of Replacing the
VerifyMatrixIdentity with the ones that have the required signature,
the changes would have been minimal .

Could you at least make clear why your changes are necessary ?

Eddy

Hi,

What I am looking at is the general improvement in performance of the benchmarks. We noticed that VerifyMatrixIdentity is one of the function hot spots in stressLinear benchmark.

On analysis, we see the code generated for TMath::Abs(m1(i,j)-m2(i,j)); is not efficient.

To check if we see better code generation and in turn better performance, I am trying to write code to directly access fElements. As you have mentioned earlier, the interpretation of the returned data pointer is different for TMatrixD, TMatrixDSym and so on. Due to this my source code changes would not work.

I was looking at changes within a function rather than have new functions defined.

Thanks,
Hrishi