From a598f67bcd2b9a1997b5ec1425215384f38aae1c Mon Sep 17 00:00:00 2001
From: Eric Cano <eric.cano@cern.ch>
Date: Thu, 8 Apr 2021 15:17:53 +0200
Subject: [PATCH] Reduced duplicate computations with temporary variables.
 (#614)

Also removed addition of transposed matrix by directly copying the needed elements.
https://github.com/cms-sw/cmssw/pull/31722#discussion_r603395679
https://github.com/cms-sw/cmssw/pull/31722#discussion_r603396371
---
 .../PixelTrackFitting/interface/BrokenLine.h  | 41 ++++++++++++-------
 1 file changed, 26 insertions(+), 15 deletions(-)

diff --git a/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h b/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h
index 86fe6a278777c..76b4bd19fc2cd 100644
--- a/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h
+++ b/RecoPixelVertexing/PixelTrackFitting/interface/BrokenLine.h
@@ -208,31 +208,42 @@ namespace brokenline {
                                                                const riemannFit::VectorNd<n>& varBeta) {
     riemannFit::MatrixNd<n> c_uMat = riemannFit::MatrixNd<n>::Zero();
     for (u_int i = 0; i < n; i++) {
+      // Pre compute differences used multiple time. Notation is t<j><k> = sTotal(i + k) - sTotal(i + k - j).
+      // We pre-compute only when possible/needed
+      // The compiler with unroll this loop with a known boundary, and no tests should
+      // be actually executed.
+      const auto t10 =  (i > 0)                   ? sTotal(i)     - sTotal(i - 1) : 0;
+      const auto t21 = ((i > 0) && (i < (n - 1))) ? sTotal(i + 1) - sTotal(i - 1) : 0;
+      const auto t11 =             (i < (n - 1))  ? sTotal(i + 1) - sTotal(i)     : 0;
+      const auto t12 =             (i < (n - 2))  ? sTotal(i + 2) - sTotal(i + 1) : 0;
+      const auto t22 =             (i < (n - 2))  ? sTotal(i + 2) - sTotal(i)     : 0;
+      const auto t11xt10 = t11 * t10;
+      const auto t12xt11 = t12 * t11;
       c_uMat(i, i) = weights(i);
       if (i > 1)
-        c_uMat(i, i) += 1. / (varBeta(i - 1) * riemannFit::sqr(sTotal(i) - sTotal(i - 1)));
+        c_uMat(i, i) += 1. / (varBeta(i - 1) * riemannFit::sqr(t10));
+
       if (i > 0 && i < n - 1)
-        c_uMat(i, i) +=
-            (1. / varBeta(i)) * riemannFit::sqr((sTotal(i + 1) - sTotal(i - 1)) /
-                                                ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1))));
+        c_uMat(i, i) += riemannFit::sqr(t21 / t11xt10) / varBeta(i);
+
       if (i < n - 2)
-        c_uMat(i, i) += 1. / (varBeta(i + 1) * riemannFit::sqr(sTotal(i + 1) - sTotal(i)));
+        c_uMat(i, i) += 1. / (varBeta(i + 1) * riemannFit::sqr(t11));
 
       if (i > 0 && i < n - 1)
-        c_uMat(i, i + 1) =
-            1. / (varBeta(i) * (sTotal(i + 1) - sTotal(i))) *
-            (-(sTotal(i + 1) - sTotal(i - 1)) / ((sTotal(i + 1) - sTotal(i)) * (sTotal(i) - sTotal(i - 1))));
-      if (i < n - 2)
-        c_uMat(i, i + 1) +=
-            1. / (varBeta(i + 1) * (sTotal(i + 1) - sTotal(i))) *
-            (-(sTotal(i + 2) - sTotal(i)) / ((sTotal(i + 2) - sTotal(i + 1)) * (sTotal(i + 1) - sTotal(i))));
+        c_uMat(i, i + 1) =  -t21 / (varBeta(i) * t11 * t11xt10);
 
       if (i < n - 2)
-        c_uMat(i, i + 2) = 1. / (varBeta(i + 1) * (sTotal(i + 2) - sTotal(i + 1)) * (sTotal(i + 1) - sTotal(i)));
+        c_uMat(i, i + 1) += -t22 / (varBeta(i + 1) * t11 * t12xt11);
 
-      c_uMat(i, i) *= 0.5;
+      if (i < n - 1) // for 2 previous steps, non diagonal element copy for this symmetric matrix
+        c_uMat(i + 1, i) = c_uMat(i, i + 1);
+
+      if (i < n - 2) {
+        c_uMat(i, i + 2) = 1. / (varBeta(i + 1) * t12xt11);
+        c_uMat(i + 2, i) = c_uMat(i, i + 2); // non-diagonal element copy.
+      }
     }
-    return c_uMat + c_uMat.transpose();
+    return c_uMat;
   }
 
   /*!