Shift of row index in TMatrixTSparse?

[berthoum]

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 [3] TMatrixT<double> Abc(Asparse)
(TMatrixT<double> &) Name: TMatrixT<double> Title: 
root [4] Abc(1,962)
(double) 0.0000000
root [5] Abc(2,962)
(double) 1.6451417e-06
root [6] RF.RowI[0]
(int) 1
root [7] RF.ColI[0]
(int) 962
root [8] RF.Aelt[0]
(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) {

instead?

Regards,
Eric Berthoumieux

---

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

---

[Eddy_Offermann]

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;
}

[Axel]

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.