arrayops.hpp

Go to the documentation of this file.

// Triumvirate: Three-Point Clustering Measurements in LSS
//
// Copyright (C) 2023 Mike S Wang & Naonori S Sugiyama [GPL-3.0-or-later]
//
// This file is part of the Triumvirate program. See the COPYRIGHT
// and LICENCE files at the top-level directory of this distribution
// for details of copyright and licensing.
//
// This program is free software: you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.

 
#ifndef TRIUMVIRATE_INCLUDE_ARRAYOPS_HPP_INCLUDED_
#define TRIUMVIRATE_INCLUDE_ARRAYOPS_HPP_INCLUDED_
 
#if defined(TRV_USE_HIP)
#include <hipfft/hipfftXt.h>
#elif defined(TRV_USE_CUDA)  // !TRV_USE_HIP && TRV_USE_CUDA
#include <cufftXt.h>
#endif                       // TRV_USE_HIP
#include <fftw3.h>
 
#include <algorithm>
#include <cfenv>
#include <cmath>
#include <cstdarg>
#include <stdexcept>
#include <vector>
 
#include "monitor.hpp"
 
#if defined(TRV_USE_HIP)
using fft_double_complex = hipfftDoubleComplex;
using fftHandle = hipfftHandle;
using libXtDesc = hipLibXtDesc;
#elif defined(TRV_USE_CUDA)  // !TRV_USE_HIP && TRV_USE_CUDA
using fft_double_complex = cufftDoubleComplex;
using fftHandle = cufftHandle;
using libXtDesc = cudaLibXtDesc;
#endif                       // TRV_USE_HIP
 
namespace trv {
 
// ***********************************************************************
// Extrapolation
// ***********************************************************************
 
namespace sys {

class ExtrapError: public std::runtime_error {
 public:
  std::string err_mesg;  

  ExtrapError(const char* fmt_string, ...);

  virtual const char* what() const noexcept;
};
 
}  // namespace trv::sys
 
namespace array {

enum ExtrapOption {
  NONE = 0,    
  ZERO = 1,    
  PAD = 2,     
  LIN = 3,     
  LOGLIN = 4,  
};

template <typename T>
bool check_isclose(
  const std::vector<T>& arr, const T val,
  T atol = 1.e-8, T rtol = 1.e-5
) {
  for (T elem : arr) {
      if (std::abs(elem - val) > atol + rtol * std::abs(val)) {return false;}
  }
  return true;
}

int check_1d_array(
  std::vector<double>& a, bool check_lin, bool check_loglin, bool check_sign
);

void extrap_lin(std::vector<double>& a, int N_ext, std::vector<double>& a_ext);

void extrap_loglin(
  std::vector<double>& a, int N_ext, std::vector<double>& a_ext
);

void extrap_pad(
  std::vector<double>& a,
  int N_ext, double c_lower, double c_upper,
  std::vector<double>& a_ext
);

void extrap2d_lin(
  std::vector< std::vector<double> >& a,
  int N_row_ext, int N_col_ext,
  std::vector< std::vector<double> >& a_ext
);

void extrap2d_loglin(
  std::vector< std::vector<double> >& a,
  int N_row_ext, int N_col_ext,
  std::vector< std::vector<double> >& a_ext
);

void extrap2d_pad(
  std::vector< std::vector<double> >& a,
  int N_row_ext, int N_col_ext,
  double c_row_lower, double c_row_upper,
  double c_col_lower, double c_col_upper,
  std::vector< std::vector<double> >& a_ext
);
 
 
// ***********************************************************************
// Sorting
// ***********************************************************************

std::vector<int> get_sorted_indices(std::vector<int> sorting_vector);
 
 
// ***********************************************************************
// Memory management
// ***********************************************************************
 
#if defined(TRV_USE_HIP) || defined(TRV_USE_CUDA)
void copy_complex_array_dtoh(
  fft_double_complex* d_arr, fftw_complex* arr, std::size_t length
);

void copy_complex_array_htod(
  fftw_complex* arr, fft_double_complex* d_arr, std::size_t length
);

void copy_complex_array_dtoh_mgpu(
  fftHandle plan, libXtDesc* libxt_desc, fftw_complex* arr
);

void copy_complex_array_htod_mgpu(
  fftHandle plan, libXtDesc* libxt_desc, fftw_complex* arr
);
#endif  // TRV_USE_HIP || TRV_USE_CUDA
}  // namespace trv::array
 
}  // namespace trv
 
#endif  // !TRIUMVIRATE_INCLUDE_ARRAYOPS_HPP_INCLUDED_