XGCa
linear_1d_interpolation.hpp
Go to the documentation of this file.
1 #ifndef LINEAR_1D_INTERPOLATION_HPP
2 #define LINEAR_1D_INTERPOLATION_HPP
3
4 #include "space_settings.hpp"
5
6
7 /* Simple 1D linear interpolation
8  * It is assumed that x is sorted in ascending order
9  * @param(in) x View of input x values
10  * @param(in) y View of input y(x) values
11  * @param(in) x1 Single x value
12  * @return T2 Single y value
13  */
14 template<typename T, typename T2>
15 KOKKOS_INLINE_FUNCTION T2 linear_1d_interpolation(const T& x, const T& y, const T2 x1){
16  const int n = x.size();
17  if (x1 >= x(0) && x1 <= x(n-1)){
18  // use bisection to find points in x closest to x1
19  int is = 0;
20  int ie = n-1;
21
22  while(ie-is > 1){
23  int ip = (is+ie)/2;
24  if (x1 > x(ip)){
25  is=ip;
26  }else{
27  ie=ip;
28  }
29  }
30
31  double slope = (y(ie)-y(is))/(x(ie)-x(is));
32  return (y(is) + slope * (x1-x(is)));
33  }else{
34  // Extrapolate if outside range
35  if (x1 < x(0)){
36  double slope = (y(1)-y(0)) / (x(1)-x(0));
37  return (y(0) + slope * (x1-x(0)));
38  }else if (x1 > x(n-1)){
39  double slope = (y(n-1)-y(n-2)) / (x(n-1)-x(n-2));
40  return (y(n-1) + slope * (x1-x(n-1)));
41  }else{
42  return 0.0;
43  // NaN?
44  }
45  }
46 }
47
48
49 /* Simple 1D linear interpolation
50  * It is assumed that x is sorted in ascending order
51  * @param(in) x View of input x values
52  * @param(in) y View of input y(x) values
53  * @param(in) x1 View of output x values
54  * @param(out) y1 View of output y(x1) values
55  */
56 template<typename T>
57 void linear_1d_interpolation(const T& x, const T& y, const T& x1, const T& y1){
58  Kokkos::parallel_for("interpol_1d", Kokkos::RangePolicy<typename T::execution_space>( 0,x1.size() ), KOKKOS_LAMBDA( const int i){
59  y1(i) = linear_1d_interpolation(x, y, x1(i));
60  });
61  Kokkos::fence();
62 }
63
64 #endif
KOKKOS_INLINE_FUNCTION T2 linear_1d_interpolation(const T &x, const T &y, const T2 x1)
Definition: linear_1d_interpolation.hpp:15
void parallel_for(const std::string name, int n_ptl, Function func, Option option, HostAoSoA aosoa_h, DeviceAoSoA aosoa_d)
Definition: streamed_parallel_for.hpp:252