Merge branch 'RemoveMathLibScalarProduct' into 'master'

Remove MathLib scalar product

See merge request ogs/ogs!3512

GeoLib/Raster.cpp

@@ -90,7 +90,9 @@ double Raster::interpolateValueAtPoint(MathLib::Point3d const& pnt) const
     // weights for bilinear interpolation
     double const xShift = std::abs((xPos - xIdx) - 0.5);
     double const yShift = std::abs((yPos - yIdx) - 0.5);
-    std::array<double,4> weight = {{ (1-xShift)*(1-yShift), xShift*(1-yShift), xShift*yShift, (1-xShift)*yShift }};
+    Eigen::Vector4d weight = {(1 - xShift) * (1 - yShift),
+                              xShift * (1 - yShift), xShift * yShift,
+                              (1 - xShift) * yShift};
 
     // neighbors to include in interpolation
     int const xShiftIdx = (xPos - xIdx >= 0.5) ? 1 : -1;
@@ -99,7 +101,7 @@ double Raster::interpolateValueAtPoint(MathLib::Point3d const& pnt) const
     std::array<int,4> const y_nb = {{ 0, 0, yShiftIdx, yShiftIdx }};
 
     // get pixel values
-    std::array<double,4>  pix_val{};
+    Eigen::Vector4d  pix_val{};
     unsigned no_data_count (0);
     for (unsigned j=0; j<4; ++j)
     {
@@ -135,12 +137,11 @@ double Raster::interpolateValueAtPoint(MathLib::Point3d const& pnt) const
             return _header.no_data;
         }
 
-        const double norm = 1.0 / (weight[0]+weight[1]+weight[2]+weight[3]);
-        std::for_each(weight.begin(), weight.end(), [&norm](double &val){val*=norm;});
+        weight /= weight.sum();
     }
 
     // new value
-    return MathLib::scalarProduct<double,4>(weight.data(), pix_val.data());
+    return weight.dot(pix_val);
     }
 
 bool Raster::isPntOnRaster(MathLib::Point3d const& pnt) const

MathLib/MathTools.h

@@ -9,60 +9,13 @@
 
 #pragma once
 
-#include <Eigen/Eigen>
 #include <cstddef>
 
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-
 namespace MathLib
 {
 template <typename T, std::size_t DIM> class TemplatePoint;
 using Point3d = MathLib::TemplatePoint<double, 3>;
 
-/**
- * standard inner product in R^N
- * \param v0 array of type T representing the vector
- * \param v1 array of type T representing the vector
- * */
-template<typename T, int N> inline
-T scalarProduct(T const * const v0, T const * const v1)
-{
-    T res (v0[0] * v1[0]);
-
-#pragma omp parallel for reduction (+:res)
-    for (int k = 1; k < N; k++)
-    {
-        res += v0[k] * v1[k];
-    }
-    return res;
-}
-
-template <> inline
-double scalarProduct<double,3>(double const * const v0, double const * const v1)
-{
-    double res (v0[0] * v1[0]);
-    for (std::size_t k(1); k < 3; k++)
-    {
-        res += v0[k] * v1[k];
-    }
-    return res;
-}
-
-template <typename T>
-inline T scalarProduct(T const* const v0, T const* const v1, int const n)
-{
-    T res (v0[0] * v1[0]);
-
-#pragma omp parallel for reduction (+:res)
-    for (int k = 1; k < n; k++)
-    {
-        res += v0[k] * v1[k];
-    }
-    return res;
-}
-
 /**
  * calcProjPntToLineAndDists computes the orthogonal projection
  * of a point p to the line described by the points a and b,