From 63613f115ee755bcde3bae115d01773e791eb147 Mon Sep 17 00:00:00 2001
From: Nicolas Morange <gartrog@melix.net>
Date: Tue, 27 Jun 2023 17:17:11 +0200
Subject: [PATCH] Fix bug in DiagProd calculation

Fixes a small bug in calculation of product of diagonal terms in LU
matrix decompositions.

A unit test to cover the former bug is also implemented.
---
 math/matrix/src/TDecompBase.cxx        |  2 +-
 math/matrix/test/CMakeLists.txt        |  1 +
 math/matrix/test/testMatrixTDecomp.cxx | 39 ++++++++++++++++++++++++++
 3 files changed, 41 insertions(+), 1 deletion(-)
 create mode 100644 math/matrix/test/testMatrixTDecomp.cxx

diff --git a/math/matrix/src/TDecompBase.cxx b/math/matrix/src/TDecompBase.cxx
index fa1d254e43912..aa6939c605789 100644
--- a/math/matrix/src/TDecompBase.cxx
+++ b/math/matrix/src/TDecompBase.cxx
@@ -228,7 +228,7 @@ void TDecompBase::DiagProd(const TVectorD &diag,Double_t tol,Double_t &d1,Double
    for (Int_t i = 0; (((i < n) && (t1 !=zero ))); i++) {
       if (TMath::Abs(diag(i)) > tol) {
          t1 *= (Double_t) diag(i);
-         while ( TMath::Abs(t1) < one) {
+         while ( TMath::Abs(t1) >= one) {
             t1 *= sixteenth;
             t2 += four;
             niter2++;
diff --git a/math/matrix/test/CMakeLists.txt b/math/matrix/test/CMakeLists.txt
index d6bc6af1d4f80..7b0739dde9d1f 100644
--- a/math/matrix/test/CMakeLists.txt
+++ b/math/matrix/test/CMakeLists.txt
@@ -5,3 +5,4 @@
 # For the list of contributors see $ROOTSYS/README/CREDITS.
 
 ROOT_ADD_GTEST(testMatrixTSparse testMatrixTSparse.cxx LIBRARIES Matrix)
+ROOT_ADD_GTEST(testMatrixTDecomp testMatrixTDecomp.cxx LIBRARIES Matrix)
diff --git a/math/matrix/test/testMatrixTDecomp.cxx b/math/matrix/test/testMatrixTDecomp.cxx
new file mode 100644
index 0000000000000..0116a7f8b825a
--- /dev/null
+++ b/math/matrix/test/testMatrixTDecomp.cxx
@@ -0,0 +1,39 @@
+// Authors: Nicolas Morange   Dec 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 <TDecompBase.h>
+
+#include <gtest/gtest.h>
+
+#include <iostream>
+
+// This is just so the can use the protected DiagProd funciton in the test.
+class TDecompDummy : public TDecompBase {
+public:
+   static void DiagProd(const TVectorD &diag, Double_t tol, Double_t &d1, Double_t &d2)
+   {
+      return TDecompBase::DiagProd(diag, tol, d1, d2);
+   }
+};
+
+// https://github.com/root-project/root/issues/13110
+TEST(testDecomp, DiagProd)
+{
+   TVectorD v(1);
+   v[0] = 1024;
+   double d1;
+   double d2;
+   TDecompDummy::DiagProd(v, 0.1, d1, d2);
+
+   // DiagProd returns the product of matrix diagonal elements in d1 and d2. d1
+   // is a mantissa and d2 an exponent for powers of 2. This is why we are
+   // using this specific formula to validate the method.
+   EXPECT_EQ(d1 * std::pow(2, d2), v[0]);
+}