# Shift of row index in TMatrixTSparse?

Hello,
I found a strange behavior by using the TMatrixTSparse class.
I created a matrix by using the code:

``````   TMatrixTSparse<double> Asparse(1,nbinsE,1,nbinsE,RowI.size(),RowI.data(),ColI.data(),Aelt.data());
``````

where `nbins=968`, `RowI` and `ColI` are `vector<int>` containing the row and column indexes of non zeroes matrix element store in the `vector<double> Aelt`

By looking the elements of the stored matrix I saw a shift of one row with respect of the input values. For instance:

``````root  TMatrixT<double> Abc(Asparse)
(TMatrixT<double> &) Name: TMatrixT<double> Title:
root  Abc(1,962)
(double) 0.0000000
root  Abc(2,962)
(double) 1.6451417e-06
root  RF.RowI
(int) 1
root  RF.ColI
(int) 962
root  RF.Aelt
(double) 1.6451417e-06
``````

By looking into the TMatrixTSparse.cxx source code (can’t post link sorry), I have the feeling there is a mistake at line 1240:

``````      if (ielem < nr && row[ielem] < irow) {
``````

Should it be

``````      if (ielem < nr && row[ielem] - this->fRowLwb< irow) {
``````

Regards,
Eric Berthoumieux

ROOT Version: 6.28.04
Platform: macosxarm64
Compiler: Apple clang 14.0.3

Hi Eric,

That is a massive bug in the initialization of a sparse matrix with a row index not starting at 0 ! This bug was able to bypass a very long list of stress tests on the matrix package in test/stressLinear.cxx

Thank you for pointing that out and debugging the problem.

Here a small reproducer code:

``````void testSparse(Int_t msize=5)
{
TMatrixDSparse m1(1,4,0,msize-1);
{
Int_t nr = 4*msize;
Int_t    *irow = new Int_t[nr];
Int_t    *icol = new Int_t[nr];
Double_t *val  = new Double_t[nr];

Int_t n = 0;
for (UInt_t i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++) {
for (UInt_t j = m1.GetColLwb(); j <= m1.GetColUpb(); j++) {
irow[n] = i;
icol[n] = j;
val[n] = TMath::Pi()*i+TMath::E()*j;
n++;
}
}
m1.SetMatrixArray(nr,irow,icol,val);
delete [] irow;
delete [] icol;
delete [] val;
}

m1.Print();

TMatrixD m2(1,4,0,msize-1);
for (UInt_t i = m2.GetRowLwb(); i <= m2.GetRowUpb(); i++)
for (UInt_t j = m2.GetColLwb(); j <= m2.GetColUpb(); j++)
m2(i,j) = TMath::Pi()*i+TMath::E()*j;
m2.Print();

std::cout << "matrices identical " << ((m1 == m2) ? "OK" : "FAILED") << std::endl;
}
``````

That’s now Incorrect initialization of TMatrixTSparse · Issue #13848 · root-project/root · GitHub . In the linked PR I have applied the patch suggested by Eric. Please leave your review comments in case I broke anything. The test is now part of ROOT’s test suite.

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.