From 323c7e48f45f3178c11d05eb6e577469aa86d348 Mon Sep 17 00:00:00 2001
From: Axel Naumann <Axel.Naumann@cern.ch>
Date: Thu, 12 Oct 2023 11:05:22 +0200
Subject: [PATCH] [matrix] Fix TMatrixTSparse::SetMatrixArray for RowLwb!=0:

Addresses https://github.com/root-project/root/issues/13848
Thanks, Eddy!
---
 math/matrix/CMakeLists.txt             |  3 ++
 math/matrix/src/TMatrixTSparse.cxx     |  4 +--
 math/matrix/test/CMakeLists.txt        |  7 ++++
 math/matrix/test/testMatrixTSparse.cxx | 46 ++++++++++++++++++++++++++
 4 files changed, 58 insertions(+), 2 deletions(-)
 create mode 100644 math/matrix/test/CMakeLists.txt
 create mode 100644 math/matrix/test/testMatrixTSparse.cxx

diff --git a/math/matrix/CMakeLists.txt b/math/matrix/CMakeLists.txt
index 7d34bcbc3a0b5..64d7e20a6d135 100644
--- 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
index 57a7b497f41a8..cf18fe2e6ea0d 100644
--- 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
new file mode 100644
index 0000000000000..d6bc6af1d4f80
--- /dev/null
+++ 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
new file mode 100644
index 0000000000000..e72dcaf803537
--- /dev/null
+++ b/math/matrix/test/testMatrixTSparse.cxx
@@ -0,0 +1,46 @@
+// 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"
+
+// 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);
+}