TMatrixD in ROOT v3 and v4

ROOT
#1

Hi,
I am trying to update a package that performs vertexing with mass/kinematic constraints. As you can imagine, it uses heavily the TMatrixD object. The package runs under ROOT v3, it compiles in v4 but the matrix manipulations fail (it appears to be an inversion but could be due to something further upstream)

Unfortunately the author is no longer in contact so I am alone in fixing things.
I have read the v4.00/08 Release Notes so know something of what has changed. These say that things are almost backward compatible except for indexing with GetElements. It appears that this member function is not used in the package.

I just wondered if anyone has done something similar, or if any ROOT experts can shed some more light on this.

Cheers
Chris

#2

Hi Chris,

I rewrote the linear algebra package . Indeed the claim is that besides
matrix storage (was column-wise, is now row-wise) , it is backward
compatible . Could you tell me a bit more how it fails :

  • wrong result ?
  • run-time error ?

If you send me the code (edmondoffermann - yahoo - com), I will
have a look at it ,

Eddy

#3

Hi Eddy, thanks for the quick response.

By fail, I mean I get the following message
if (ifail) cout << “W Matrix inversion failed!” << endl;
in the code included below.

I see manipulations such as
new_winv[i][j]
could this be the problem?

I can also send the verbose output if needed,

Cheers
Chris

TMatrixD winv(TMatrixD(TMatrixD::kTransposed,a)
,TMatrixD::kMult
,TMatrixD(v
,TMatrixD::kMult
,a));

Int_t nrows = winv.GetNrows();
HepSymMatrix new_winv(nrows);
for (Int_t i=0;i<nrows;i++){
for (Int_t j=0;j<nrows;j++){
   new_winv[i][j] = winv(i,j);
}
}
Int_t ifail;
new_winv.invert(ifail);
if (ifail) cout << "W Matrix inversion failed!" << endl;
for (Int_t i=0;i<nrows;i++){
for (Int_t j=0;j<nrows;j++){
   winv(i,j) = new_winv[i][j];
}
}
derr = 1.0;
if (ifail) derr=0.0;
//winv.Invert(&derr);
w = winv;

#ifdef VERBOSE
cout << "the inverted W matrix is: " << endl;
for (Int_t i=0;i<nconstraints;i++){
for (Int_t j=0;j<nconstraints;j++){
cout << w(i,j) << " ";
}
cout << endl;
}
cout << endl;
#endif

#4

Hi Chris,

The signature of the invert is :

TMatrixD &TMatrixD::Invert(Double_t *det)

You are using :

Int_t ifail;
new_winv.invert(ifail);

So I am not surprised that the “cerr << …” is triggered . I do
not understand that that was not happening in root version 3 !!

Since I see some CLHEP stuff floating around I guess that someone
changed :

void HepMatrix::invert(int &ierr) NOTE THE int !!!

to a ROOT expression and forgot to change the argument .

So make the change and run it again

Eddy

#5

Hi,

Int_t ifail;
new_winv.invert(ifail);

In the code, new_winv is a HepSymMatrix so I think it has the correct syntax. For some reason a number of matrix inversions are done using the HepSymMatrix class, ie a TMatrixD is converted to a HepSymMatrix, inverted and then converted back again. I have no idea why this is the case. I changed the code to just use TMatrixD but keep getting the same problem. It seems the matrix that I am trying to invert is just junk (det zero).
The cause is still under investigation, the one thing I know it that it worked in v3!

Cheers
Chris

#6

Hi Chris,

Indeed, It is a CLHEP matrix so the syntax is right . ROOT
is only used then for somple things like transpose, matrix
multiplications .

It would be interesting to print a matrix content when the inversion
fails in Root 4 but not in 3 .

Eddy

#7

Hi Eddy,
On further investigation, I see the same behaviour in both v3 and v4 of ROOT now I have built both in gcc. The older version was built with kcc and I think this maybe the difference.

Cheers
Chris