KfUtils.h

Go to the documentation of this file.

/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Sergey Gorbunov, Sergei Zharko [committer] */
 
 
#pragma once  // include this header only once per compilation unit
 
#include "KfDefs.h"
#include "KfSimd.h"
 
#include <array>
#include <sstream>
 
namespace cbm::algo::kf::utils
{
 
  template<typename T>
  using scaltype = typename std::conditional<std::is_same<T, fvec>::value, fscal, T>::type;
 
  template<typename T>
  using masktype = typename std::conditional<std::is_same<T, fvec>::value, fmask, bool>::type;
 
  inline fvec iif(const fmask& m, const fvec& t, const fvec& f) { return Vc::iif(m, t, f); }
  inline fvec fabs(const fvec& v) { return Vc::abs(v); }
 
  template<typename T>
  inline T iif(bool b, T t, T f)
  {
    return b ? t : f;
  }
 
  template<typename T>
  inline T fabs(const T& v)
  {
    return std::fabs(v);
  }
 
  template<typename T>
  inline bool isFull(const T& m)
  {
    if constexpr (std::is_same_v<T, fmask>) {
      return m.isFull();
    }
    else {
      return m;
    }
  }
 
  template<typename T>
  inline T max(const T& a, const T& b)
  {
    if constexpr (std::is_same_v<T, fvec>) {
      return Vc::max(a, b);
    }
    else {
      return std::max(a, b);
    }
  }
 
#if defined(__GNUC__) && __GNUC__ - 0 >= 3
  __attribute__((always_inline)) inline bool IsLikely(bool b) { return __builtin_expect(!!(b), 1); }
  __attribute__((always_inline)) inline bool IsUnlikely(bool b) { return __builtin_expect(!!(b), 0); }
#else
  inline bool IsLikely(bool b) { return b; }
  inline bool IsUnlikely(bool b) { return b; }
#endif
 
  template<typename T>
  inline bool IsUndefined(const T& val)
  {
    if constexpr (std::is_same_v<T, float>) {
      return std::isnan(val);
    }
    else if constexpr (std::is_same_v<T, double>) {
      return std::isnan(val);
    }
    else if constexpr (std::is_same_v<T, fscal>) {
      return std::isnan(val);
    }
    else if constexpr (std::is_same_v<T, fvec>) {
      return isnan(val).isNotEmpty();
    }
    else {
      return val == cbm::algo::kf::defs::Undef<T>;
    }
  }
 
  template<typename T>
  inline bool IsFinite(const T& val)
  {
    if constexpr (std::is_same_v<T, float>) {
      return std::isfinite(val);
    }
    else if constexpr (std::is_same_v<T, double>) {
      return std::isfinite(val);
    }
    else if constexpr (std::is_same_v<T, fscal>) {
      return std::isfinite(val);
    }
    else if constexpr (std::is_same_v<T, fvec>) {
      return isfinite(val).isFull();
    }
    else {
      return val != cbm::algo::kf::defs::Undef<T>;
    }
  }
 
  template<typename DataT>
  void PrintSIMDmsg(std::stringstream& msg, DataT& v);
 
  template<typename DataT>
  inline void CheckSimdVectorEquality(DataT /*v*/, const char* /*name*/)
  {
  }
 
  template<>
  void CheckSimdVectorEquality(fvec v, const char* name);
 
  template<typename TdataA, typename TdataB, bool TDoAllA, bool TDoAllB>
  class VecCopySpec {
   public:
    [[gnu::always_inline]] static void CopyEntries(TdataA& a, int ia, const TdataB& b, int ib);
  };
 
  template<typename TdataA, typename TdataB>
  class VecCopySpec<TdataA, TdataB, true, true> {
   public:
    [[gnu::always_inline]] static void CopyEntries(TdataA& a, int, const TdataB& b, int) { a = b; }
  };
 
  template<typename TdataA, typename TdataB>
  class VecCopySpec<TdataA, TdataB, true, false> {
   public:
    [[gnu::always_inline]] static void CopyEntries(TdataA& a, int, const TdataB& b, int ib) { a = b[ib]; }
  };
 
  template<typename TdataA, typename TdataB>
  class VecCopySpec<TdataA, TdataB, false, true> {
   public:
    [[gnu::always_inline]] static void CopyEntries(TdataA& a, int ia, const TdataB& b, int) { a[ia] = b; }
  };
 
  template<typename TdataA, typename TdataB>
  class VecCopySpec<TdataA, TdataB, false, false> {
   public:
    [[gnu::always_inline]] static void CopyEntries(TdataA& a, int ia, const TdataB& b, int ib) { a[ia] = b[ib]; }
  };
 
 
  template<typename TdataA, typename TdataB, bool TDoAllA, bool TDoAllB>
  class VecCopy {
   public:
    static void CopyEntries(TdataA& a, int ia, const TdataB& b, int ib)
    {
      constexpr bool allA{TDoAllA || !std::is_same<TdataA, fvec>::value};
      constexpr bool allB{TDoAllB || !std::is_same<TdataB, fvec>::value};
      VecCopySpec<TdataA, TdataB, allA, allB>::CopyEntries(a, ia, b, ib);
    }
  };
 
  struct EnumClassHash {
    template<typename T>
    int operator()(T t) const
    {
      return static_cast<int>(t);
    }
  };
}  // namespace cbm::algo::kf::utils
 
namespace cbm::algo::kf::utils::simd
{
  template<typename DataT, typename DataOut>
  inline DataOut Cast(const DataT& val)
  {
    return static_cast<DataOut>(val);
  }
 
  template<>
  inline float Cast(const fvec& val)
  {
    return val[0];
  }
 
  template<>
  inline double Cast(const fvec& val)
  {
    return static_cast<double>(val[0]);
  }
 
  template<typename DataT, typename DataOut>
  inline DataOut Cast(const DataT& val, size_t /*i*/)
  {
    return static_cast<DataOut>(val);
  }
 
  template<typename DataT, typename DataOut, size_t N>
  inline std::array<DataOut, N> Cast(const std::array<DataT, N>& arr, size_t i)
  {
    std::array<DataOut, N> res;
    for (size_t iEl = 0; iEl < arr.size(); ++iEl) {
      res[iEl] = Cast<DataT, DataOut>(arr[iEl], i);
    }
    return res;
  }
 
  template<typename DataT, typename DataOut, size_t N>
  inline std::array<DataOut, N> Cast(const std::array<DataT, N>& arr)
  {
    std::array<DataOut, N> res;
    for (size_t iEl = 0; iEl < arr.size(); ++iEl) {
      res[iEl] = Cast<DataT, DataOut>(arr[iEl]);
    }
    return res;
  }
 
  template<>
  inline float Cast(const fvec& val, size_t i)
  {
    return val[i];
  }
 
  template<>
  inline double Cast(const fvec& val, size_t i)
  {
    return static_cast<double>(val[i]);
  }
 
  template<typename DataT>
  inline size_t Size()
  {
    return 1;
  }
 
  template<>
  inline size_t Size<fvec>()
  {
    return fvec::size();
  }
 
  template<typename DataT>
  inline DataT One()
  {
    return static_cast<DataT>(1);
  }
 
  template<>
  inline fvec One<fvec>()
  {
    return fvec::One();
  }
 
  template<typename DataT, typename DataIn>
  inline void SetEntry(DataT& out, DataIn in, size_t /*i*/)
  {
    out = static_cast<DataT>(in);
  }
 
  template<typename DataIn>
  inline void SetEntry(fvec& out, DataIn in, size_t i)
  {
    out[i] = in;
  }
}  // namespace cbm::algo::kf::utils::simd
 
 
// TODO: SZh 06.06.2024:
//       There is nothing to do with SIMD in the funcitons below. Thus, we should move them into KfMath.h
namespace cbm::algo::kf::utils::math
{
 
  void CholeskyFactorization(const double a[], const int n, int nn, double u[], int* nullty, int* ifault);
 
  void SymInv(const double a[], const int n, double c[], double w[], int* nullty, int* ifault);
 
}  // namespace cbm::algo::kf::utils::math