fftlog.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_FFTLOG_HPP_INCLUDED_
#define TRIUMVIRATE_INCLUDE_FFTLOG_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 <cmath>
#include <complex>
#include <cstring>
#include <vector>
 
#include "maths.hpp"
#include "arrayops.hpp"
 
#if defined(TRV_USE_HIP)
using fft_double_complex = hipfftDoubleComplex;
using fftHandle = hipfftHandle;
#elif defined(TRV_USE_CUDA)  // !TRV_USE_HIP && TRV_USE_CUDA
using fft_double_complex = cufftDoubleComplex;
using fftHandle = cufftHandle;
#endif                       // TRV_USE_HIP
 
namespace trva = trv::array;
 
namespace trv {
 
namespace maths {

class HankelTransform {
 public:
  double order;         
  double bias;          
  int nsamp = 0;        
  int nsamp_trans = 0;  
  double logres = 0.;   
  double pivot = 1.;    

  std::vector<double> pre_sampts;

  std::vector<double> post_sampts;

  HankelTransform(double mu, double q, bool threaded = true);

  ~HankelTransform();

  void initialise(
    std::vector<double> sample_pts, double kr_c,
    bool lowring = true,
    trva::ExtrapOption extrap = trva::ExtrapOption::NONE,
    double extrap_exp = 2.
  );

  void initialise(
    std::vector<double> sample_pts, double kr_c,
    bool lowring = true,
    int extrap = 0,
    double extrap_exp = 2.
  );

  void biased_transform(std::complex<double>* a, std::complex<double>* b);

  std::vector< std::complex<double> > compute_kernel_coeff();

  double calc_lowring_pivot(double delta, double kr_c = 1.);
 
 private:
  trva::ExtrapOption extrap = trva::ExtrapOption::NONE;

  int n_ext = 0;

  std::vector<double> pre_sampts_extrap;

  std::vector<double> post_sampts_extrap;

  std::vector< std::complex<double> > kernel;

  fftw_plan pre_plan{};
  fftw_plan post_plan{};
#if defined(TRV_USE_HIP) || defined(TRV_USE_CUDA)
  fftHandle plan_gpu{};
#endif  // TRV_USE_HIP || TRV_USE_CUDA

  fftw_complex* pre_buffer = nullptr;
  fftw_complex* post_buffer = nullptr;
#if defined(TRV_USE_HIP) || defined(TRV_USE_CUDA)
  fft_double_complex* d_buffer = nullptr;
#endif  // TRV_USE_HIP || TRV_USE_CUDA
 
#ifdef _CUDA_STREAM
  cudaStream_t custream{};
#endif  // _CUDA_STREAM

  bool plan_init = false;

  bool threaded = true;

  void reset();
};
 
class SphericalBesselTransform: public HankelTransform {
 public:
  int degree;  

  SphericalBesselTransform(int ell, int n, bool threaded = true);

  void initialise(
    std::vector<double> sample_pts, double kr_c,
    bool lowring = true,
    trva::ExtrapOption extrap = trva::ExtrapOption::NONE,
    double extrap_exp = 2.
  );

  void initialise(
    std::vector<double> sample_pts, double kr_c,
    bool lowring = true,
    int extrap = 0,
    double extrap_exp = 2.
  );

  void biased_transform(
    std::vector< std::complex<double> >& a,
    std::vector< std::complex<double> >& b
  );

  void transform_cosmological_multipole(
    int dir,
    std::vector< std::complex<double> >& pre_samples,
    std::vector< std::complex<double> >& post_samples
  );
};
 
}  // namespace trv::maths
 
}  // namespace trv
 
#endif  // !TRIUMVIRATE_INCLUDE_FFTLOG_HPP_INCLUDED_