[matrix] Fix TMatrixTSparse::SetMatrixArray for RowLwb!=0:

Addresses root-project#13848 Thanks, Eddy!
Axel-Naumann · Oct 27, 2023 · 9ff6706 · 9ff6706
1 parent f5040c3
commit 9ff6706
Show file tree

Hide file tree

Showing 4 changed files with 60 additions and 2 deletions.
diff --git a/math/matrix/CMakeLists.txt b/math/matrix/CMakeLists.txt
@@ -81,3 +81,6 @@ ROOT_STANDARD_LIBRARY_PACKAGE(Matrix
  DICTIONARY_OPTIONS
    -writeEmptyRootPCM
 )
+
+ROOT_ADD_TEST_SUBDIRECTORY(test)
+
diff --git a/math/matrix/src/TMatrixTSparse.cxx b/math/matrix/src/TMatrixTSparse.cxx
@@ -1185,7 +1185,7 @@ TMatrixTBase<Element> &TMatrixTSparse<Element>::SetMatrixArray(Int_t nr,Int_t *r
    R__ASSERT(row[irowmin] >= this->fRowLwb && row[irowmax] <= this->fRowLwb+this->fNrows-1);
    R__ASSERT(col[icolmin] >= this->fColLwb && col[icolmax] <= this->fColLwb+this->fNcols-1);
 
-   if (row[irowmin] < this->fRowLwb || row[irowmax] > this->fRowLwb+this->fNrows-1) {
+   if (row[irowmin] <  this->fRowLwb|| row[irowmax] > this->fRowLwb+this->fNrows-1) {
       Error("SetMatrixArray","Inconsistency between row index and its range");
       if (row[irowmin] < this->fRowLwb) {
          Info("SetMatrixArray","row index lower bound adjusted to %d",row[irowmin]);
@@ -1237,7 +1237,7 @@ TMatrixTBase<Element> &TMatrixTSparse<Element>::SetMatrixArray(Int_t nr,Int_t *r
    Int_t ielem = 0;
    nr_nonzeros = 0;
    for (Int_t irow = 1; irow < this->fNrows+1; irow++) {
-      if (ielem < nr && row[ielem] < irow) {
+      if (ielem < nr && row[ielem] - this->fRowLwb < irow) {
          while (ielem < nr) {
             if (data[ielem] != 0.0) {
                fColIndex[nr_nonzeros] = col[ielem]-this->fColLwb;

diff --git a/math/matrix/test/CMakeLists.txt b/math/matrix/test/CMakeLists.txt
@@ -0,0 +1,7 @@
+# Copyright (C) 1995-2023, Rene Brun and Fons Rademakers.
+# All rights reserved.
+#
+# For the licensing terms see $ROOTSYS/LICENSE.
+# For the list of contributors see $ROOTSYS/README/CREDITS.
+
+ROOT_ADD_GTEST(testMatrixTSparse testMatrixTSparse.cxx LIBRARIES Matrix)
diff --git a/math/matrix/test/testMatrixTSparse.cxx b/math/matrix/test/testMatrixTSparse.cxx
@@ -0,0 +1,48 @@
+// Authors: Eddy Offermann   Oct 2023
+
+/*************************************************************************
+ * Copyright (C) 1995-2023, Rene Brun and Fons Rademakers.               *
+ * All rights reserved.                                                  *
+ *                                                                       *
+ * For the licensing terms see $ROOTSYS/LICENSE.                         *
+ * For the list of contributors see $ROOTSYS/README/CREDITS.             *
+ *************************************************************************/
+
+#include "TMatrixD.h"
+#include "TMatrixDSparse.h"
+#include "TMath.h"
+
+#include "gtest/gtest.h"
+
+#include <array>
+
+// https://github.com/root-project/root/issues/13848
+TEST(testSparse, LwbInit)
+{
+  constexpr int msize = 5;
+  TMatrixDSparse m1(1, 4, 0, msize - 1);
+  {
+    constexpr int nr = 4 * msize;
+    std::array<int, nr> irow;
+    std::array<int, nr> icol;
+    std::array<double, nr> val;
+
+    Int_t n = 0;
+    for (int i = m1.GetRowLwb(); i <= m1.GetRowUpb(); i++) {
+      for (int 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.data(), icol.data(), val.data());
+  }
+
+  TMatrixD m2(1, 4, 0, msize - 1);
+  for (int i = m2.GetRowLwb(); i <= m2.GetRowUpb(); i++)
+    for (int j = m2.GetColLwb(); j <= m2.GetColUpb(); j++)
+      m2(i,j) = TMath::Pi() * i + TMath::E() * j;
+
+  EXPECT_EQ(m1, m2);
+}