Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_cubic_splines.F90
Go to the documentation of this file.
1 !**************************************************************
2 ! Copyright INRIA
3 ! Authors :
4 ! CALVI project team
5 !
6 ! This code SeLaLib (for Semi-Lagrangian-Library)
7 ! is a parallel library for simulating the plasma turbulence
8 ! in a tokamak.
9 !
10 ! This software is governed by the CeCILL-B license
11 ! under French law and abiding by the rules of distribution
12 ! of free software. You can use, modify and redistribute
13 ! the software under the terms of the CeCILL-B license as
14 ! circulated by CEA, CNRS and INRIA at the following URL
15 ! "http://www.cecill.info".
16 !**************************************************************
17 
31 
33 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34 #include "sll_assert.h"
35 #include "sll_errors.h"
36 #include "sll_memory.h"
37 #include "sll_working_precision.h"
38 
40  sll_p_hermite, &
41  sll_p_periodic
42 
43  use sll_m_tridiagonal, only: &
46 
47  implicit none
48 
49  public :: &
50  ! 1D
60  sll_f_cubic_spline_1d_get_x1_delta, &
61  ! 2D
71  sll_f_cubic_spline_2d_get_x1_min, &
72  sll_f_cubic_spline_2d_get_x1_max, &
73  sll_f_cubic_spline_2d_get_x1_delta, &
74  sll_f_cubic_spline_2d_get_x2_min, &
75  sll_f_cubic_spline_2d_get_x2_max, &
76  sll_f_cubic_spline_2d_get_x2_delta
77 
78  private
79 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
80 
81  ! --- compile-time constants to avoid run-time division
82  sll_real64, parameter :: inv_6 = 1._f64/6._f64
83 
89  sll_int32, private :: n_points
90  sll_real64, private :: delta
91  sll_real64, private :: rdelta
92  sll_real64, private :: xmin
93  sll_real64, private :: xmax
94  sll_int32, private :: bc_type
97  sll_real64, dimension(:), pointer, private :: d => null()
99  sll_real64, dimension(:), pointer, private :: coeffs => null()
101  sll_real64, private :: slope_l
103  sll_real64, private :: slope_r
105  logical, private :: compute_slope_l
107  logical, private :: compute_slope_r
111  logical, private :: use_fast_algorithm
113  sll_real64, dimension(:), pointer, private :: a => null()
115  sll_real64, dimension(:), pointer, private :: cts => null()
117  sll_int32, dimension(:), pointer, private :: ipiv => null()
119  sll_real64, dimension(:), pointer, private :: f_aux => null()
120  end type sll_t_cubic_spline_1d
121 
128  sll_int32, private :: num_pts_x1
129  sll_int32, private :: num_pts_x2
130  sll_real64, private :: x1_delta
131  sll_real64, private :: x1_rdelta
132  sll_real64, private :: x2_delta
133  sll_real64, private :: x2_rdelta
134  sll_real64, private :: x1_min
135  sll_real64, private :: x1_max
136  sll_real64, private :: x2_min
137  sll_real64, private :: x2_max
138  sll_int32, private :: x1_bc_type
139  sll_int32, private :: x2_bc_type
140  ! if data is not used, it should be deleted make a decision...
141  sll_real64, pointer, private :: data(:, :) => null()
142  sll_real64, pointer, private :: d1(:) => null()
145  sll_real64, pointer, private :: d2(:) => null()
146  sll_real64, pointer, private :: coeffs(:, :) => null()
147  sll_real64, pointer, private :: x1_min_slopes(:) => null()
148  sll_real64, pointer, private :: x1_max_slopes(:) => null()
149  sll_real64, pointer, private :: x2_min_slopes(:) => null()
150  sll_real64, pointer, private :: x2_max_slopes(:) => null()
151  sll_real64, pointer, private :: x1_min_slopes_coeffs(:) => null()
152  sll_real64, pointer, private :: x1_max_slopes_coeffs(:) => null()
153  sll_real64, pointer, private :: x2_min_slopes_coeffs(:) => null()
154  sll_real64, pointer, private :: x2_max_slopes_coeffs(:) => null()
155  logical, private :: compute_slopes_x1_min
156  logical, private :: compute_slopes_x1_max
157  logical, private :: compute_slopes_x2_min
158  logical, private :: compute_slopes_x2_max
159  end type sll_t_cubic_spline_2d
160 
161  ! The following was written after declaring the contents of the fundamental
162  ! spline types as private. This may be an overkill. If this ever gives
163  ! problems, it is suggested to change the corresponding slots to 'public'.
164 
165 #ifndef DOXYGEN_SHOULD_SKIP_THIS
166 
167  ! Some useful macros that should probably be put in a different file to
168  ! make them more widely accessible. Here we use them to compute default
169  ! values for the slopes.
170  ! Actually the proper way to factor this is with a function that takes
171  ! an array segment of a certain length and returns the slope. Such
172  ! function(s) could be made to work in 1D or 2D by passing the right
173  ! array segment. Keep this in mind...
174 
175  ! (-25/12, 4, -3, 4/3, -1/4) stencil
176 #define FORWARD_FD_5PT( f, r_delta ) \
177  r_delta*(-(25.0_f64/12.0_f64)*f(1) + 4.0_f64*f(2) - 3.0_f64*f(3) + (4.0_f64/3.0_f64)*f(4) - 0.25_f64*f(5))
178 
179  ! (1/4, -4/3, 3, -4, 25/12) stencil
180 #define BACKWARD_FD_5PT( f, r_delta, np ) \
181  r_delta*(0.25_f64*f(np - 4) - (4.0_f64/3.0_f64)*f(np - 3) + 3.0_f64*f(np - 2) - 4.0_f64*f(np - 1) + (25.0_f64/12.0_f64)*f(np))
182 
183  ! (-3/2, 2, -1/2 ) stencil
184 #define FORWARD_FD_3PT( f, r_delta ) \
185  r_delta*(-1.5_f64*f(1) + 2.0_f64*f(2) - 0.5_f64*f(3))
186 
187  ! ( 1/2, -2, 3/2 ) stencil
188 #define BACKWARD_FD_3PT( f, r_delta, np ) \
189  r_delta*(0.5_f64*f(np - 2) - 2.0_f64*f(np - 1) + 1.5_f64*f(np))
190 
191 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
192 
193 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
194 contains
195 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
196 
197 #ifndef DOXYGEN_SHOULD_SKIP_THIS
198 
199 #define MAKE_GET_SLOT_FUNCTION( fname, datatype, slot, ret_type ) \
200  function fname(spline_obj) result(val); \
201  type(datatype) :: spline_obj; \
202  ret_type :: val; \
203  val = spline_obj%slot; \
204  end function fname
205 
206  make_get_slot_function(sll_f_cubic_spline_1d_get_x1_min, sll_t_cubic_spline_1d, xmin, sll_real64)
207  make_get_slot_function(sll_f_cubic_spline_1d_get_x1_max, sll_t_cubic_spline_1d, xmax, sll_real64)
208  make_get_slot_function(sll_f_cubic_spline_1d_get_x1_delta, sll_t_cubic_spline_1d, delta, sll_real64)
209 
210  make_get_slot_function(sll_f_cubic_spline_2d_get_x1_min, sll_t_cubic_spline_2d, x1_min, sll_real64)
211  make_get_slot_function(sll_f_cubic_spline_2d_get_x1_max, sll_t_cubic_spline_2d, x1_max, sll_real64)
212  make_get_slot_function(sll_f_cubic_spline_2d_get_x2_min, sll_t_cubic_spline_2d, x2_min, sll_real64)
213  make_get_slot_function(sll_f_cubic_spline_2d_get_x2_max, sll_t_cubic_spline_2d, x2_max, sll_real64)
214  make_get_slot_function(sll_f_cubic_spline_2d_get_x1_delta, sll_t_cubic_spline_2d, x1_delta, sll_real64)
215  make_get_slot_function(sll_f_cubic_spline_2d_get_x2_delta, sll_t_cubic_spline_2d, x2_delta, sll_real64)
216 
217  ! The following implementation embodies the algorithm described in
218  ! Eric Sonnendrucker's "A possibly faster algorithm for cubic splines on
219  ! a uniform grid" (unpublished).
220 
221  ! The array of spline coefficients has NP+3 elements. The extra elements
222  ! at the ends (i.e.: 0, NP+1, NP+2) store coefficients whose values are
223  ! determined by the type of boundary condition used. This is invisible
224  ! to the user, who should not be concerned with this implementation detail.
225 
226  ! The following parameter determines the problem size at which the alternative
227  ! spline algorithm is used.
228 
229 #define NUM_TERMS 27
230 
231 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
232 
249  subroutine sll_s_cubic_spline_1d_init(self, num_points, xmin, xmax, bc_type, sl, sr, fast_algorithm)
250  type(sll_t_cubic_spline_1d) :: self
251  sll_int32, intent(in) :: num_points
252  sll_real64, intent(in) :: xmin
253  sll_real64, intent(in) :: xmax
254  sll_int32, intent(in) :: bc_type
255  sll_real64, intent(in), optional :: sl
256  sll_real64, intent(in), optional :: sr
257  logical, intent(in), optional :: fast_algorithm
258  sll_int32 :: ierr
259  sll_int32 :: i
260 
261  self%n_points = num_points
262  self%xmin = xmin
263  self%xmax = xmax
264  self%delta = (xmax - xmin)/real((num_points - 1), f64)
265  self%rdelta = 1.0_f64/self%delta
266  self%bc_type = bc_type
267  if (num_points .lt. num_terms) then
268  self%use_fast_algorithm = .false.
269  else
270  self%use_fast_algorithm = .true.
271  if (present(fast_algorithm)) then
272  self%use_fast_algorithm = fast_algorithm
273  end if
274  end if
275  if (xmin .gt. xmax) then
276  print *, 'ERROR, sll_s_cubic_spline_1d_init: xmin is greater than xmax, ', &
277  'this would cause all sorts of errors.'
278  stop
279  end if
280  ! Some more general error checking depending on the type of boundary
281  ! condition requested.
282  select case (bc_type)
283  case (sll_p_periodic)
284  if (present(sl) .or. present(sr)) then
285 ! print *, 'self(): it is not allowed to specify the ',&
286 ! 'end ifin the case of periodic boundary conditions. ', &
287 ! 'Exiting program...'
288 ! STOP 'self'
289  sll_warning('self', 'values of sl and sr are not taken into account')
290  else
291  ! Assign some value, but this value should never be used in the
292  ! periodic case anyway.
293  self%slope_L = 0.0_f64
294  self%slope_R = 0.0_f64
295  end if
296  if (self%use_fast_algorithm .eqv. .false.) then
297  sll_allocate(self%a(3*(num_points - 1)), ierr)
298  sll_allocate(self%cts(7*(num_points - 1)), ierr)
299  sll_allocate(self%ipiv(num_points - 1), ierr)
300  self%f_aux => null() ! not needed in periodic case
301  ! Initialize and factorize the tridiagonal system. See detailed
302  ! comment below regarding the structure of this matrix.
303  self%a(1) = 1.0_f64/6.0_f64
304  self%a(2) = 4.0_f64/6.0_f64
305  self%a(3) = 1.0_f64/6.0_f64
306  self%a(3*num_points - 5) = 1.0_f64/6.0_f64
307  self%a(3*num_points - 4) = 4.0_f64/6.0_f64
308  self%a(3*num_points - 3) = 1.0_f64/6.0_f64
309  do i = 1, num_points - 3
310  self%a(3*i + 1) = 1.0_f64/6.0_f64
311  self%a(3*i + 2) = 4.0_f64/6.0_f64
312  self%a(3*i + 3) = 1.0_f64/6.0_f64
313  end do
315  self%a, &
316  num_points - 1, &
317  self%cts, &
318  self%ipiv)
319  else
320  self%a => null()
321  self%cts => null()
322  self%ipiv => null()
323  self%f_aux => null()
324  end if
325  case (sll_p_hermite)
326  if (present(sl)) then
327  self%slope_L = sl
328  self%compute_slope_L = .false.
329  else
330  self%slope_L = 0.0_f64 ! just a filler value
331  self%compute_slope_L = .true.
332  end if
333  if (present(sr)) then
334  self%slope_R = sr
335  self%compute_slope_R = .false.
336  else
337  self%slope_R = 0.0_f64 ! just a filler value
338  self%compute_slope_R = .true.
339  end if
340  if (self%use_fast_algorithm .eqv. .false.) then
341  sll_allocate(self%a(3*(num_points + 2)), ierr)
342  sll_allocate(self%cts(7*(num_points + 2)), ierr)
343  sll_allocate(self%ipiv(num_points + 2), ierr)
344  sll_allocate(self%f_aux(num_points + 2), ierr)
345  ! Initialize and factorize the tridiagonal system. See detailed
346  ! comment below regarding the structure of this matrix. Matrix 'a'
347  ! is (np+2)X(np+2)
348  self%a(1) = 0.0_f64
349  self%a(2) = 4.0_f64/6.0_f64
350  self%a(3) = 2.0_f64/6.0_f64
351  self%a(3*(num_points + 2) - 2) = 2.0_f64/6.0_f64
352  self%a(3*(num_points + 2) - 1) = 4.0_f64/6.0_f64
353  self%a(3*(num_points + 2)) = 0.0_f64
354  do i = 1, num_points
355  self%a(3*i + 1) = 1.0_f64/6.0_f64
356  self%a(3*i + 2) = 4.0_f64/6.0_f64
357  self%a(3*i + 3) = 1.0_f64/6.0_f64
358  end do
360  self%a, &
361  num_points + 2, &
362  self%cts, &
363  self%ipiv)
364  else
365  self%a => null()
366  self%cts => null()
367  self%ipiv => null()
368  self%f_aux => null()
369  end if
370  case default
371  print *, 'ERROR: sll_s_cubic_spline_1d_init(): not recognized boundary condition'
372  stop
373  end select
374  if (self%use_fast_algorithm .eqv. .true.) then
375  sll_allocate(self%d(num_points), ierr)
376  else
377  self%d => null()
378  end if
379  ! note how the indexing of the coefficients array includes the end-
380  ! points 0, num_points, num_points+1, num_points+2. These are meant to
381  ! store the boundary condition-specific data. The 'periodic' BC does
382  ! not use the num_points+2 point.
383  sll_allocate(self%coeffs(0:num_points + 2), ierr)
384  end subroutine sll_s_cubic_spline_1d_init
385 
386  ! Fast spline algorithm description:
387  !
388  ! - data: the array whose data must be fit with the cubic spline.
389  ! - np: (number of points, length of the data array that must be fit with
390  ! the spline.
391  ! - bc_type: an integer flag describing the type of boundary conditions
392  ! desired.
393  ! - spline_obj: the spline object to be initialized.
394  ! This version assumes that the data are uniformly spaced.
395  !
396  ! The idea behind this spline implementation is the factorization of the
397  ! array:
398  !
399  ! + 4 1 1 +
400  ! | 1 4 1 |
401  ! | . . . |
402  ! | . . . |
403  ! | . . . |
404  ! | 1 4 1 |
405  ! + 1 1 4 +
406  !
407  ! in the form:
408  !
409  ! A = L*L^t
410  !
411  ! where:
412  !
413  ! + a b +
414  ! | b a |
415  ! | . . |
416  ! L = | . . | (zeros not shown)
417  ! | . . |
418  ! | b a |
419  ! + b a +
420  !
421  ! This factorization is achieved for the values:
422  !
423  ! a² = (2+sqrt(3))/6, b² = (2-sqrt(3))/6.
424  !
425  ! The spline coefficients C are thus the solution of:
426  !
427  ! L*L^t*C = F.
428  !
429  ! Hence, we solve first for D in:
430  !
431  ! L*D = F
432  !
433  ! This solution is achieved in a two-step process. First we compute the
434  ! first term d1, which is given by:
435  !
436  ! d1 =
437  !
438  ! 1 b b 2 i b i
439  !---(f(x1) - ---f(xN) + (---)*f(xN-1) + ... + (-1)(---) (f(xN-i+1)-bd(N-i)))
440  ! a a a a
441  !
442  ! The series converges (since a > b) and we can approximate the series
443  ! by a partial sum:
444  !
445  ! 1 b
446  ! d1 = ---(f(x1) + SUM(i=1,M)(-1)^i(---)^i*f(xN-i+1))
447  ! a a
448  !
449  ! The rest of the terms can be found by:
450  !
451  ! d(i) = 1/a*(f(x i) - b*d(i-1))
452  !
453  ! Once D is known, the same procedure can be used to compute C in
454  !
455  ! L^t*C = D
456  !
457  ! The same process is carried out. First, the last coefficient is
458  ! calculated by:
459  !
460  ! c(N) = 1/a*(d(N) + SUM(i=1,M) (-1)^i*(b/a)^i*d(i))
461  !
462  ! And the rest of the terms, starting with C(N-1) and working backwards:
463  !
464  ! c(i) = 1/a*(d(i) - b*c(i+1))
465  !
466  ! The algorithm above is not implemented whenever the number of points is
467  ! smaller than 28. In such cases we fall back to a more straightforward but
468  ! also more costly implementation using a standard tridiagonal system
469  ! solver.
470  !
471  ! In the periodic case, we are solving the problem A*C=F, where
472  !
473  ! + 4 1 1 + 1
474  ! | 1 4 1 |
475  ! 1 | . . . | .
476  ! A = --- | . . . | .
477  ! 6 | . . . | .
478  ! | 1 4 1 |
479  ! + 1 1 4 + num_points
480  !
481  ! In the Hermite BC case, the problem is modified as follows:
482  !
483  ! + 4 2 + + c_0 + + 6*f_0+2*delta*f_0´ + 0
484  ! | 1 4 1 | | | | 6*f_1 | 1
485  ! | . . . | | . | | . | .
486  ! A =| . . . |*| . | = | . | .
487  ! | . . . | | . | | . | .
488  ! | 1 4 1 | | | | 6*f_np | np
489  ! + 2 4 + +c_(np+1)+ + 6*f_(np+1)+2*delta*f_(np+1)'+ np+1
490 
499  sll_real64, dimension(:), intent(in) :: f
500  sll_int32 :: bc_type
501  type(sll_t_cubic_spline_1d) :: spline
502  ! Note that this function does no error checking and basically
503  ! outsources this task to the functions it is wrapping around.
504  ! This is so because those functions can be used independently
505  ! (if the user wants to avoid the overhead of calling this
506  ! wrapper function), so in any case, the error checking of
507  ! the arguments will be carried out at least once.
508  bc_type = spline%bc_type;
509  select case (bc_type)
510  case (sll_p_periodic)
511  call compute_spline_1d_periodic(f, spline)
512  case (sll_p_hermite)
513  call compute_spline_1d_hermite(f, spline)
514  case default
515  print *, 'ERROR: compute_spline_1D(): not recognized boundary condition'
516  stop
517  end select
519 
520  ! The following auxiliary functions:
521  ! compute_spline_1D_periodic_aux()
522  ! compute_spline_1D_hermite_aux()
523  ! are the fundamental building blocks. These are meant to do the work
524  ! needed to compute the splines. Other functions are essentially
525  ! wrappers around these. Clients of these routines are responsible for
526  ! all error-checking.
527  ! f: data for which the cubic spline fit is desired.
528  ! num_pts: number of points that compose the data
529  ! d: scratch array, size num_pts, needed to compute the coefficients.
530  ! coeffs: output of the computation, size 0:num_pts+2
531  subroutine compute_spline_1d_periodic_aux(f, num_pts, d, coeffs)
532  sll_real64, dimension(:), pointer :: f
533  sll_int32, intent(in) :: num_pts
534  sll_real64, dimension(:), pointer :: d
535  sll_real64, dimension(:), pointer :: coeffs
536  sll_real64, parameter :: a = sqrt((2.0_f64 + sqrt(3.0_f64))/6.0_f64)
537  sll_real64, parameter :: r_a = 1.0_f64/a
538  sll_real64, parameter :: b = sqrt((2.0_f64 - sqrt(3.0_f64))/6.0_f64)
539  sll_real64, parameter :: b_a = b/a
540  sll_real64 :: coeff_tmp
541  sll_real64 :: d1
542  sll_int32 :: i
543  sll_int32 :: np
544  sll_assert(size(f) .ge. num_pts - 1)
545  sll_assert(size(d) .ge. num_pts)
546  sll_assert(size(coeffs) .ge. num_pts)
547  sll_assert(num_pts .gt. num_terms)
548 
549  np = num_pts
550  ! Compute d(1):
551  d1 = f(1)
552  coeff_tmp = 1.0_f64
553  do i = 0, num_terms - 1 ! if NUM_TERMS == 0, only f(nc) is considered.
554 ! do i = 1, NUM_TERMS
555  coeff_tmp = coeff_tmp*(-b_a)
556  d1 = d1 + coeff_tmp*f(np - 1 - i)
557  end do
558  ! Fill the d array with the intermediate result
559  d(1) = d1*r_a
560  do i = 2, np - 1
561  d(i) = r_a*(f(i) - b*d(i - 1))
562  end do
563  ! Compute the coefficients. Start with first term
564  d1 = d(np - 1)
565  coeff_tmp = 1.0_f64
566  do i = 1, num_terms
567  coeff_tmp = coeff_tmp*(-b_a)
568  d1 = d1 + coeff_tmp*d(i)
569  end do
570 ! coeffs(np-1) = d1*r_a
571  coeffs(np) = d1*r_a
572  ! rest of the coefficients:
573  do i = np - 2, 1, -1
574 ! coeffs(i) = r_a*(d(i) - b*coeffs(i+1))
575  coeffs(i + 1) = r_a*(d(i) - b*coeffs(i + 2))
576  end do
577  coeffs(1) = coeffs(np)
578  coeffs(np + 1) = coeffs(2)
579  coeffs(np + 2) = coeffs(3)
580  coeffs(np + 3) = coeffs(4)
581  end subroutine compute_spline_1d_periodic_aux
582 
584  f, &
585  num_pts, &
586  d, &
587  slope_l, &
588  slope_r, &
589  delta, &
590  coeffs)
591 
592  sll_real64, dimension(:), pointer :: f
593  sll_int32, intent(in) :: num_pts
594  sll_real64, dimension(:), pointer :: d
595  sll_real64, intent(in) :: slope_l
596  sll_real64, intent(in) :: slope_r
597  sll_real64, intent(in) :: delta
598  sll_real64, dimension(:), pointer :: coeffs
599  sll_int32 :: i
600  sll_int32 :: np
601  sll_real64, parameter :: a = sqrt((2.0_f64 + sqrt(3.0_f64))/6.0_f64)
602  sll_real64, parameter :: r_a = 1.0_f64/a
603  sll_real64, parameter :: b = sqrt((2.0_f64 - sqrt(3.0_f64))/6.0_f64)
604  sll_real64, parameter :: b_a = b/a
605  sll_real64, parameter :: ralpha = sqrt(6.0_f64/sqrt(3.0_f64))
606  sll_real64 :: coeff_tmp
607  sll_real64 :: d1
608  sll_real64 :: f1 ! to store modified value of f(1)
609  sll_real64 :: fnp ! for modified value of f(np)
610  np = num_pts
611  ! For Hermitian boundary conditions with non-zero slope, we can use the
612  ! same algorithm than for the zero-slope case, with the difference that
613  ! we tweak the values of the source term, f, right at the endpoints.
614  !
615  ! This can be derived from the formula for S'(x_N), which yields
616  !
617  ! S'(x_N) = 1/(2*h)*[-C_(N-1) + C_(N+1)]
618  !
619  ! For a slope = 0, this reduces to C_(N+1) = C_(N-1). For slope != 0, we
620  ! have an extra constant term (2*h*S') that we can absorb with the
621  ! source term.
622  f1 = f(1)
623  fnp = f(np) - delta*slope_r/3.0_f64
624 
625  ! Compute d(1):
626  d1 = f1
627  coeff_tmp = 1.0_f64
628 
629  ! Since we want to consider the case in which there is a given slope
630  ! at point 1, we assume that all points to the left of point 1 are
631  ! displaced linearly with an offset given by a line with the given slope.
632  ! This requires the subtraction of the 2*slope*delta*(i-1) term...
633  do i = 2, num_terms ! if NUM_TERMS == 2, only f(2) is considered.
634  coeff_tmp = coeff_tmp*(-b_a)
635  d1 = d1 + coeff_tmp*(f(i) - 2.0_f64*slope_l*delta*real(i - 1, f64))
636  end do
637  ! Fill the d array with the intermediate result
638  d(1) = d1*r_a
639  do i = 2, np - 1
640  d(i) = r_a*(f(i) - b*d(i - 1))
641  end do
642  d(np) = ralpha*(0.5_f64*fnp - b*d(np - 1))
643  ! Compute the coefficients. Start with first term
644 ! coeffs(np) = ralpha*d(np)
645  coeffs(np + 1) = ralpha*d(np)
646  do i = np - 1, 1, -1
647  coeffs(i + 1) = r_a*(d(i) - b*coeffs(i + 2))
648  end do
649  coeffs(1) = coeffs(3) - 2.0_f64*delta*slope_l
650  coeffs(np + 2) = coeffs(np) + 2.0_f64*delta*slope_r
651  coeffs(np + 3) = 0.0_f64 !coeffs(np-2) ! not used
652  end subroutine compute_spline_1d_hermite_aux
653 
654  subroutine compute_spline_1d_periodic(f, spline)
655  sll_real64, dimension(:), intent(in), target :: f ! data to be fit
656  type(sll_t_cubic_spline_1d) :: spline
657  sll_real64, dimension(:), pointer :: coeffs
658  sll_int32 :: np
659  sll_real64, dimension(:), pointer :: d
660  sll_real64, dimension(:), pointer :: fp
661 
662  if (.not. (size(f) .ge. spline%n_points - 1)) then
663  ! FIXME: THROW ERROR
664  print *, 'ERROR: compute_spline_1D_periodic(): '
665  write (*, '(a, i8, a, i8)') 'spline object needs data of size >= ', &
666  spline%n_points - 1, ' . Passed size: ', size(f)
667  stop
668  end if
669  if (spline%use_fast_algorithm .eqv. .false.) then
670  np = spline%n_points - 1 ! remove periodic point for this algorithm
672  spline%cts, &
673  spline%ipiv, &
674  f, &
675  np, &
676  spline%coeffs(1:np))
677  ! and set the periodic BC by setting the coefficient values
678  spline%coeffs(0) = spline%coeffs(np)
679  spline%coeffs(np + 1) = spline%coeffs(1)
680  spline%coeffs(np + 2) = spline%coeffs(2)
681  spline%coeffs(np + 3) = spline%coeffs(3)
682  else
683  np = spline%n_points
684  fp => f
685  d => spline%d
686  coeffs => spline%coeffs(0:np + 2)
687  ! Remember that now coeffs(1) refers to spline%coeffs(0)
688  call compute_spline_1d_periodic_aux(fp, np, d, coeffs)
689  end if
690  end subroutine compute_spline_1d_periodic
691 
692  subroutine compute_spline_1d_hermite(f, spline)
693  sll_real64, dimension(:), intent(in), target :: f ! data to be fit
694  type(sll_t_cubic_spline_1d) :: spline
695  sll_real64, dimension(:), pointer :: coeffs
696  sll_int32 :: np
697  sll_real64, dimension(:), pointer :: fp
698  sll_real64, dimension(:), pointer :: d
699  sll_real64 :: slope_l
700  sll_real64 :: slope_r
701  sll_real64 :: delta
702  sll_real64 :: r_delta ! reciprocal
703  !sll_int32 :: i
704 
705  if (.not. (size(f) .ge. spline%n_points)) then
706  ! FIXME: THROW ERROR
707  print *, 'ERROR: compute_spline_1D_hermite(): '
708  write (*, '(a, i8, a, i8)') 'spline object needs data of size >= ', &
709  spline%n_points, ' . Passed size: ', size(f)
710  stop
711  end if
712  fp => f
713  np = spline%n_points
714  d => spline%d
715  coeffs => spline%coeffs(0:)
716  delta = spline%delta
717  r_delta = 1.0_f64/delta
718 
719  ! Estimate numerically the values of the slopes based on the given
720  ! values of 'f' if the user did not provide values initially.
721  if (spline%compute_slope_L .eqv. .true.) then
722  slope_l = forward_fd_5pt(f, r_delta)
723  else
724  slope_l = spline%slope_L
725  end if
726  if (spline%compute_slope_R .eqv. .true.) then
727  slope_r = backward_fd_5pt(f, r_delta, np)
728  else
729  slope_r = spline%slope_R
730  end if
731 
732  if (.not. spline%use_fast_algorithm) then
733  ! load the source term.
734  spline%f_aux(2:np + 1) = fp(1:np)
735  ! set the end-points to reflect the slope information
736  spline%f_aux(1) = fp(1) + (1.0_f64/3.0_f64)*delta*slope_l
737  spline%f_aux(np + 2) = fp(np) - (1.0_f64/3.0_f64)*delta*slope_r
739  spline%cts, &
740  spline%ipiv, &
741  spline%f_aux, &
742  np + 2, &
743  spline%coeffs)
744  else
746  fp, np, d, slope_l, slope_r, delta, coeffs)
747  end if
748  end subroutine compute_spline_1d_hermite
749 
750  ! Here we use the cubic B-spline centered at node 'i', supported on
751  ! four cells, two on each side of node 'i':
752  !
753  ! S_(i-1) S_i S_(i+1) S_(i+2)
754  ! + + + + +
755  ! | | | | |
756  !
757  ! *
758  ! * *
759  ! * *
760  ! * *
761  ! * *
762  ! * *
763  ! * *
764  ! --------*--------+--------+--------+--------*-----
765  ! i-2 i-1 i i+1 i+2
766  !
767  ! From left to right, the curves that define the above shape piecewise are:
768  !
769  ! S_(i-1)(x) defined in [x_(i-2), x_(i-1)]:
770  !
771  ! = - 1/(6h^3)*(x_(i-2)-x)^3
772  !
773  ! S_i(x) defined in [x_(i-1), x_i]:
774  !
775  ! = 1/(6h^3)*[3*(x_(i-1)-x)^3 + 3*h*(x_(i-1)-x)^2 - 3*h^2*(x_(i-1)-x) + h^3]
776  !
777  ! S_(i+1)(x) defined in [x_i, x_(i+1)]:
778  !
779  ! = 1/(6h^3)*[3*(x-x_(i+1))^3 + 3*h*(x-x_(i+1))^2 - 3*h^2*(x-x_(i+1)) + h^3]
780  !
781  ! S_(i+2)(x) defined in [x_(i+1), x_(i+2)]:
782  !
783  ! = - 1/(6h^3)*(x - x_(i+2))^3
784  !
785  ! Thus any given point will have a value defined by the 4 contributions of
786  ! the splines that are supported in the same interval as the point.
787  ! In this implementation, we choose a normalized representation, in which
788  ! h = 1 and a given point is represented by the number of the cell it lives
789  ! in and the local offset within this cell. The offset is also normalized to
790  ! unity.
791  !
792  ! The image of a given point 'x', under a function modeled by the splines
793  ! would be:
794  !
795  ! S(x) = C_(i-1)*S_(i+2) + C_i*S_(i+1) + C_(i+1)*S_i + C_(i+2)*S_(i-1)
796  !
797  ! In normalized form (h=1), this can be written:
798  !
799  ! S(x) =
800  !
801  ! 1/6*C_(i-1)*(1-dx)^3 + 1/6*C_i*[-3*(1-dx)^3 + 3*(1-dx)^2 +3*(1-dx) + 1] +
802  !
803  ! 1/6*C_(i+1)*[-3*dx^3 + 3*dx^2 + 3*dx + 1] + 1/6*C_(i+2)*dx^3
804  !
805  ! The above is the actual interpolation. Here we attempt to minimize the
806  ! number of operations by carrying out a few factorizations. Hence, what is
807  ! implemented may be something more akin to:
808  !
809  ! 1/6*[(1-dx)*[(1-dx)*[(1-dx)*[C_(i-1) - 3*C_i] + 3*C_i] + 3*C_i] + C_i +
810  !
811  ! dx*[dx*[dx*[C_(i+2) - 3*C_(i+1)] + 3*C_(i+1)] + 3*C_(i+1)] + C_(i+1)]
812  !
813  ! As a check, when the point has a zero offset (dx = 0), the above formula
814  ! reduces correctly to:
815  !
816  ! 1/6*[C_(i-1) + 4*C_i + C_(i+1)]
817  !
818  ! Similarly, if the offset is unity (dx = 1), the above formula reduces to
819  !
820  ! 1/6*[C_i + 4*C_(i+1) + C_(i+2)]
821  !
822  ! as needed.
823 
824  ! interpolate_value_aux() is meant to be used as a private function. It
825  ! aims at providing a reusable code to carry out interpolations in the
826  ! other module functions (like higher dimensions). This is why this
827  ! function does no argument checking. Arguments should be checked by
828  ! the caller.
829  function interpolate_value_aux(x, xmin, rh, coeffs)
830  sll_real64 :: interpolate_value_aux
831  sll_real64, intent(in) :: x
832  sll_real64, intent(in) :: xmin
833  sll_real64, dimension(:), pointer :: coeffs
834  sll_real64, intent(in) :: rh ! reciprocal of cell spacing
835  sll_int32 :: cell
836  sll_real64 :: dx
837  sll_real64 :: cdx ! 1-dx
838  sll_real64 :: t0 ! temp/scratch variables ...
839  sll_real64 :: t1
840  sll_real64 :: t2
841  sll_real64 :: t3
842  sll_real64 :: t4
843  sll_real64 :: cim1 ! C_(i-1)
844  sll_real64 :: ci ! C_i
845  sll_real64 :: cip1 ! C_(i+1)
846  sll_real64 :: cip2 ! C_(i+2)
847  ! find the cell and offset for x
848  t0 = (x - xmin)*rh
849  cell = int(t0) + 1
850  dx = t0 - real(cell - 1, f64)
851  cdx = 1.0_f64 - dx
852  cim1 = coeffs(cell)
853  ci = coeffs(cell + 1)
854  cip1 = coeffs(cell + 2)
855  cip2 = coeffs(cell + 3)
856  ! print *, 'intepolate_value_aux(): coefficients:'
857  ! print *, cim1, ci, cip1, cip2, cdx, dx
858  t1 = 3.0_f64*ci
859  t3 = 3.0_f64*cip1
860  t2 = cdx*(cdx*(cdx*(cim1 - t1) + t1) + t1) + ci
861  t4 = dx*(dx*(dx*(cip2 - t3) + t3) + t3) + cip1
862  ! print *, 't2 and t4: ', t2, t4
863  interpolate_value_aux = (1.0_f64/6.0_f64)*(t2 + t4)
864  !print *, interpolate_value_aux
865  end function interpolate_value_aux
866 
874  function sll_f_cubic_spline_1d_eval(x, spline)
875  sll_real64 :: sll_f_cubic_spline_1d_eval
876  intrinsic :: associated, int, real
877  sll_real64, intent(in) :: x
878  type(sll_t_cubic_spline_1d) :: spline
879  sll_real64, dimension(:), pointer :: coeffs
880  sll_real64 :: xmin
881  sll_real64 :: rh ! reciprocal of cell spacing
882 
883  ! We set these as assertions since we want the flexibility of turning
884  ! them off.
885  sll_assert((x .ge. spline%xmin) .and. (x .le. spline%xmax))
886  xmin = spline%xmin
887  rh = spline%rdelta
888  coeffs => spline%coeffs(0:spline%n_points + 2)
890  end function sll_f_cubic_spline_1d_eval
891 
903  subroutine sll_s_cubic_spline_1d_eval_array(a_in, a_out, n, spline)
904  intrinsic :: associated, int, real
905  sll_int32, intent(in) :: n
906  sll_real64, dimension(1:n), intent(in) :: a_in
907  sll_real64, dimension(1:n), intent(out) :: a_out
908  type(sll_t_cubic_spline_1d) :: spline
909  sll_real64, dimension(:), pointer :: coeffs
910  sll_real64 :: rh ! reciprocal of cell spacing
911  sll_int32 :: cell
912  sll_real64 :: dx
913  sll_real64 :: cdx ! 1-dx
914  sll_real64 :: t0 ! temp/scratch variables ...
915  sll_real64 :: t1
916  sll_real64 :: t2
917  sll_real64 :: t3
918  sll_real64 :: t4
919  sll_real64 :: cim1 ! C_(i-1)
920  sll_real64 :: ci ! C_i
921  sll_real64 :: cip1 ! C_(i+1)
922  sll_real64 :: cip2 ! C_(i+2)
923  sll_int32 :: num_cells
924  sll_real64 :: x
925  sll_int32 :: i
926  ! FIXME: arg checks here
927  num_cells = spline%n_points - 1
928  rh = spline%rdelta
929  coeffs => spline%coeffs
930  ! find the cell and offset for x
931  do i = 1, n
932  x = a_in(i)
933  if (.not. ((x .ge. spline%xmin) .and. (x .le. spline%xmax))) then
934  print *, 'splines', x, spline%xmin, spline%xmax
935  end if
936  !print*, 'splines', x, spline%xmin, spline%xmax
937  sll_assert((x .ge. spline%xmin) .and. (x .le. spline%xmax))
938  t0 = (x - spline%xmin)*rh
939  cell = int(t0) + 1
940  dx = t0 - real(cell - 1, f64)
941  cdx = 1.0_f64 - dx
942  ! write (*,'(a,i15, a, f20.12)') 'cell = ', cell, ', dx = ', dx
943  cim1 = coeffs(cell - 1)
944  ci = coeffs(cell)
945  cip1 = coeffs(cell + 1)
946  cip2 = coeffs(cell + 2)
947  t1 = 3.0_f64*ci
948  t3 = 3.0_f64*cip1
949  t2 = cdx*(cdx*(cdx*(cim1 - t1) + t1) + t1) + ci
950  t4 = dx*(dx*(dx*(cip2 - t3) + t3) + t3) + cip1
951  a_out(i) = (1.0_f64/6.0_f64)*(t2 + t4)
952  !print*,'sll_s_cubic_spline_1d_eval_array', i, a_out(i)
953  end do
954  end subroutine sll_s_cubic_spline_1d_eval_array
955 
956  ! FIXME: The following function is not in the unit test.
957 
971  subroutine interpolate_pointer_values(ptr_in, ptr_out, n, spline)
972  intrinsic :: associated, int, real
973  sll_int32, intent(in) :: n
974  sll_real64, dimension(:), pointer :: ptr_in
975  sll_real64, dimension(:), pointer :: ptr_out
976  type(sll_t_cubic_spline_1d) :: spline
977  sll_real64, dimension(:), pointer :: coeffs
978  sll_real64 :: rh ! reciprocal of cell spacing
979  sll_int32 :: cell
980  sll_real64 :: dx
981  sll_real64 :: cdx ! 1-dx
982  sll_real64 :: t0 ! temp/scratch variables ...
983  sll_real64 :: t1
984  sll_real64 :: t2
985  sll_real64 :: t3
986  sll_real64 :: t4
987  sll_real64 :: cim1 ! C_(i-1)
988  sll_real64 :: ci ! C_i
989  sll_real64 :: cip1 ! C_(i+1)
990  sll_real64 :: cip2 ! C_(i+2)
991  sll_int32 :: num_cells
992  sll_real64 :: x
993  sll_int32 :: i
994  sll_assert(associated(ptr_in))
995  sll_assert(associated(ptr_out))
996  ! FIXME: arg checks here
997  num_cells = spline%n_points - 1
998  rh = spline%rdelta
999  coeffs => spline%coeffs
1000  ! find the cell and offset for x
1001  do i = 1, n
1002  x = ptr_in(i)
1003  !print*, 'splines', x, spline%xmin, spline%xmax
1004  sll_assert((x .ge. spline%xmin) .and. (x .le. spline%xmax))
1005  t0 = (x - spline%xmin)*rh
1006  cell = int(t0) + 1
1007  dx = t0 - real(cell - 1, f64)
1008  cdx = 1.0_f64 - dx
1009  ! write (*,'(a,i8, a, f20.12)') 'cell = ', cell, ', dx = ', dx
1010  cim1 = coeffs(cell - 1)
1011  ci = coeffs(cell)
1012  cip1 = coeffs(cell + 1)
1013  cip2 = coeffs(cell + 2)
1014  t1 = 3.0_f64*ci
1015  t3 = 3.0_f64*cip1
1016  t2 = cdx*(cdx*(cdx*(cim1 - t1) + t1) + t1) + ci
1017  t4 = dx*(dx*(dx*(cip2 - t3) + t3) + t3) + cip1
1018  ptr_out(i) = (1.0_f64/6.0_f64)*(t2 + t4)
1019  end do
1020  end subroutine interpolate_pointer_values
1021 
1022  ! interpolate_derivative_aux() is a private function aimed at abstracting
1023  ! away the capability of computing the derivative at a point, given the
1024  ! array of cubic spline coefficients.
1025  function interpolate_derivative_aux(x, xmin, rh, coeffs)
1026  sll_real64 :: interpolate_derivative_aux
1027  intrinsic :: int, real
1028  sll_real64, intent(in) :: x
1029  sll_real64, intent(in) :: xmin
1030  sll_real64, intent(in) :: rh ! reciprocal of cell spacing
1031  sll_real64, dimension(:), pointer :: coeffs
1032  sll_int32 :: cell
1033  sll_real64 :: dx
1034  sll_real64 :: t0 ! temp/scratch variables ...
1035  sll_real64 :: t1
1036  sll_real64 :: t2
1037  sll_real64 :: t3
1038  sll_real64 :: cim1 ! C_(i-1)
1039  sll_real64 :: ci ! C_i
1040  sll_real64 :: cip1 ! C_(i+1)
1041  sll_real64 :: cip2 ! C_(i+2)
1042 
1043  ! find the cell and offset for x
1044  t0 = (x - xmin)*rh
1045  cell = int(t0) + 1
1046  dx = t0 - real(cell - 1, f64)
1047  ! write (*,'(a,i8, a, f20.12)') 'cell = ', cell, ', dx = ', dx
1048  cim1 = coeffs(cell)
1049  ci = coeffs(cell + 1)
1050  cip1 = coeffs(cell + 2)
1051  cip2 = coeffs(cell + 3)
1052  t1 = 2.0_f64*(cim1 - 2.0_f64*ci + cip1)
1053  t2 = -cim1 + 3.0_f64*(ci - cip1) + cip2
1054  t3 = cip1 - cim1
1055  interpolate_derivative_aux = 0.5_f64*rh*(dx*(t1 + dx*t2) + t3)
1056  end function interpolate_derivative_aux
1057 
1066  sll_real64 :: sll_f_cubic_spline_1d_eval_deriv
1067  intrinsic :: associated
1068  sll_real64, intent(in) :: x
1069  sll_real64, dimension(:), pointer :: coeffs
1070  type(sll_t_cubic_spline_1d) :: spline
1071 
1072  ! We set these as assertions since we want the flexibility of turning
1073  ! them off.
1074  sll_assert((x .ge. spline%xmin) .and. (x .le. spline%xmax))
1075  coeffs => spline%coeffs(0:spline%n_points + 2)
1077  x, &
1078  spline%xmin, &
1079  spline%rdelta, &
1080  coeffs)
1082 
1095  array_in, &
1096  array_out, &
1097  num_pts, &
1098  spline)
1099 
1100  intrinsic :: associated
1101  sll_real64, dimension(:), intent(in) :: array_in
1102  sll_int32, intent(in) :: num_pts
1103  sll_real64, dimension(:), intent(out) :: array_out
1104  type(sll_t_cubic_spline_1d) :: spline
1105  sll_real64, dimension(:), pointer :: coeffs
1106  sll_int32 :: i
1107 
1108  sll_assert(num_pts .le. size(array_in))
1109 
1110  coeffs => spline%coeffs(0:spline%n_points + 2)
1111  do i = 1, num_pts
1112  sll_assert((array_in(i) .ge. spline%xmin) .and. (array_in(i) .le. spline%xmax))
1113  array_out(i) = interpolate_derivative_aux( &
1114  array_in(i), spline%xmin, spline%rdelta, coeffs)
1115  end do
1117 
1118  ! FIXME: The following subroutine is not in the unit test
1129  ptr_in, &
1130  ptr_out, &
1131  num_pts, &
1132  spline)
1133 
1134  intrinsic :: associated
1135  sll_real64, dimension(:), pointer :: ptr_in
1136  sll_int32, intent(in) :: num_pts
1137  sll_real64, dimension(:), pointer :: ptr_out
1138  type(sll_t_cubic_spline_1d) :: spline
1139  sll_real64, dimension(:), pointer :: coeffs
1140  sll_int32 :: i
1141 
1142  sll_assert(num_pts .le. size(ptr_in))
1143  sll_assert(associated(ptr_in))
1144  sll_assert(associated(ptr_out))
1145  coeffs => spline%coeffs(0:spline%n_points + 2)
1146 
1147  do i = 1, num_pts
1148  sll_assert((ptr_in(i) .ge. spline%xmin) .and. (ptr_in(i) .le. spline%xmax))
1149  ptr_out(i) = interpolate_derivative_aux( &
1150  ptr_in(i), spline%xmin, spline%rdelta, coeffs)
1151  end do
1152  end subroutine interpolate_pointer_derivatives
1153 
1155  subroutine sll_s_cubic_spline_1d_free(spline)
1156  type(sll_t_cubic_spline_1d) :: spline
1157  sll_int32 :: ierr
1158 
1159  if (spline%use_fast_algorithm .eqv. .true.) then
1160  sll_deallocate(spline%d, ierr)
1161  end if
1162  sll_deallocate(spline%coeffs, ierr)
1163 ! spline%data => null()
1164  if (spline%use_fast_algorithm .eqv. .false.) then
1165  sll_deallocate(spline%a, ierr)
1166  sll_deallocate(spline%cts, ierr)
1167  sll_deallocate(spline%ipiv, ierr)
1168  if (spline%bc_type == sll_p_hermite) then
1169  sll_deallocate(spline%f_aux, ierr)
1170  end if
1171  end if
1172  end subroutine sll_s_cubic_spline_1d_free
1173 
1174  !-----------------------------------------------------------------------
1175  !
1176  ! Functions and subroutines for the 2D spline.
1177  !
1178  !----------------------------------------------------------------------
1179 
1230  num_pts_x1, &
1231  num_pts_x2, &
1232  x1_min, &
1233  x1_max, &
1234  x2_min, &
1235  x2_max, &
1236  x1_bc_type, &
1237  x2_bc_type, &
1238  const_slope_x1_min, &
1239  const_slope_x1_max, &
1240  const_slope_x2_min, &
1241  const_slope_x2_max, &
1242  x1_min_slopes, &
1243  x1_max_slopes, &
1244  x2_min_slopes, &
1245  x2_max_slopes) result(this)
1246 
1247  type(sll_t_cubic_spline_2d), pointer :: this
1248  sll_int32, intent(in) :: num_pts_x1
1249  sll_int32, intent(in) :: num_pts_x2
1250  sll_real64, intent(in) :: x1_min
1251  sll_real64, intent(in) :: x1_max
1252  sll_real64, intent(in) :: x2_min
1253  sll_real64, intent(in) :: x2_max
1254  sll_int32, intent(in) :: x1_bc_type
1255  sll_int32, intent(in) :: x2_bc_type
1256  sll_real64, intent(in), optional :: const_slope_x1_min
1257  sll_real64, intent(in), optional :: const_slope_x1_max
1258  sll_real64, intent(in), optional :: const_slope_x2_min
1259  sll_real64, intent(in), optional :: const_slope_x2_max
1260  sll_real64, intent(in), dimension(:), optional :: x1_min_slopes
1261  sll_real64, intent(in), dimension(:), optional :: x1_max_slopes
1262  sll_real64, intent(in), dimension(:), optional :: x2_min_slopes
1263  sll_real64, intent(in), dimension(:), optional :: x2_max_slopes
1264 
1265  allocate (this)
1267  this, &
1268  num_pts_x1, &
1269  num_pts_x2, &
1270  x1_min, &
1271  x1_max, &
1272  x2_min, &
1273  x2_max, &
1274  x1_bc_type, &
1275  x2_bc_type, &
1276  const_slope_x1_min, &
1277  const_slope_x1_max, &
1278  const_slope_x2_min, &
1279  const_slope_x2_max, &
1280  x1_min_slopes, &
1281  x1_max_slopes, &
1282  x2_min_slopes, &
1283  x2_max_slopes)
1284 
1285  end function sll_f_new_cubic_spline_2d
1286 
1288  this, &
1289  num_pts_x1, &
1290  num_pts_x2, &
1291  x1_min, &
1292  x1_max, &
1293  x2_min, &
1294  x2_max, &
1295  x1_bc_type, &
1296  x2_bc_type, &
1297  const_slope_x1_min, &
1298  const_slope_x1_max, &
1299  const_slope_x2_min, &
1300  const_slope_x2_max, &
1301  x1_min_slopes, &
1302  x1_max_slopes, &
1303  x2_min_slopes, &
1304  x2_max_slopes)
1305 
1306  type(sll_t_cubic_spline_2d) :: this
1307  sll_int32, intent(in) :: num_pts_x1
1308  sll_int32, intent(in) :: num_pts_x2
1309  sll_real64, intent(in) :: x1_min
1310  sll_real64, intent(in) :: x1_max
1311  sll_real64, intent(in) :: x2_min
1312  sll_real64, intent(in) :: x2_max
1313  sll_int32, intent(in) :: x1_bc_type
1314  sll_int32, intent(in) :: x2_bc_type
1315  sll_real64, intent(in), optional :: const_slope_x1_min
1316  sll_real64, intent(in), optional :: const_slope_x1_max
1317  sll_real64, intent(in), optional :: const_slope_x2_min
1318  sll_real64, intent(in), optional :: const_slope_x2_max
1319  sll_real64, intent(in), dimension(:), optional :: x1_min_slopes
1320  sll_real64, intent(in), dimension(:), optional :: x1_max_slopes
1321  sll_real64, intent(in), dimension(:), optional :: x2_min_slopes
1322  sll_real64, intent(in), dimension(:), optional :: x2_max_slopes
1323 
1324  sll_int32 :: bc_selector
1325  sll_int32 :: ierr
1326 
1327  this%num_pts_x1 = num_pts_x1
1328  this%num_pts_x2 = num_pts_x2
1329  this%x1_min = x1_min
1330  this%x1_max = x1_max
1331  this%x2_min = x2_min
1332  this%x2_max = x2_max
1333  this%x1_delta = (x1_max - x1_min)/real((num_pts_x1 - 1), f64)
1334  this%x2_delta = (x2_max - x2_min)/real((num_pts_x2 - 1), f64)
1335  this%x1_rdelta = 1.0_f64/this%x1_delta
1336  this%x2_rdelta = 1.0_f64/this%x2_delta
1337  this%x1_bc_type = x1_bc_type
1338  this%x2_bc_type = x2_bc_type
1339 
1340  if ((num_pts_x1 .le. num_terms) .or. (num_pts_x2 .le. num_terms)) then
1341  sll_error( "sll_s_cubic_spline_2d_init"," Because of the algorithm used, this function is meant to be used with arrays that are at least of size = 28")
1342  end if
1343  if ((x1_min .gt. x1_max) .or. (x2_min .gt. x2_max)) then
1344  sll_error( "sll_s_cubic_spline_2d_init"," one of the xmin is greater than the corresponding xmax, this would cause all sorts of errors")
1345  end if
1346 
1347  ! Check that slope arrays are of the right size. Consider making this
1348  ! something more permanent than an assertion. Does this compile if the
1349  ! assertions are turned off???
1350  if (present(x1_min_slopes)) then
1351  sll_assert(size(x1_min_slopes) .ge. num_pts_x2)
1352  end if
1353  if (present(x1_max_slopes)) then
1354  sll_assert(size(x1_max_slopes) .ge. num_pts_x2)
1355  end if
1356  if (present(x2_min_slopes)) then
1357  sll_assert(size(x2_min_slopes) .ge. num_pts_x1)
1358  end if
1359  if (present(x2_max_slopes)) then
1360  sll_assert(size(x2_max_slopes) .ge. num_pts_x1)
1361  end if
1362 
1363  sll_allocate(this%d1(num_pts_x1), ierr)
1364  sll_allocate(this%d2(num_pts_x2), ierr)
1365 
1366  ! Treat the bc_selector variable essentially like a bit field, to
1367  ! accumulate the information on the different boundary conditions
1368  ! given. This scheme allows to add more types of boundary conditions
1369  ! if necessary.
1370  bc_selector = 0
1371  if (x1_bc_type .eq. sll_p_periodic) then
1372  bc_selector = bc_selector + 1
1373  end if
1374  if (x1_bc_type .eq. sll_p_hermite) then
1375  bc_selector = bc_selector + 2
1376  end if
1377  if (x2_bc_type .eq. sll_p_periodic) then
1378  bc_selector = bc_selector + 4
1379  end if
1380  if (x2_bc_type .eq. sll_p_hermite) then
1381  bc_selector = bc_selector + 8
1382  end if
1383 
1384  select case (bc_selector)
1385  case (5)
1386  ! both boundary condition types are periodic
1387  if ( &
1388  present(x1_min_slopes) .or. present(x1_max_slopes) .or. &
1389  present(x2_min_slopes) .or. present(x2_max_slopes) .or. &
1390  present(const_slope_x1_min) .or. present(const_slope_x1_max) .or. &
1391  present(const_slope_x2_min) .or. present(const_slope_x2_max)) then
1392 
1393  sll_warning('sll_s_cubic_spline_2d_init', 'values of slopes are not taken into account as we are in periodic-periodic')
1394  end if
1395  case (6)
1396  ! Hermite condition in X1 and periodic in X2
1397  if ( &
1398  present(x2_min_slopes) .or. present(x2_max_slopes) .or. &
1399  present(const_slope_x2_min) .or. present(const_slope_x2_max)) then
1400  sll_warning('sll_s_cubic_spline_2d_init', 'values of slopes in x2 are not taken into account as we are periodic in x2')
1401 
1402  end if
1403  if (present(const_slope_x1_min) .and. present(x1_min_slopes)) then
1404  sll_error("sll_s_cubic_spline_2d_init"," hermite-periodic-case, it is not allowed to specify simultaneously a constant value for the slopes at x1_min and an array-specified set of slopes")
1405  end if
1406  if (present(const_slope_x1_max) .and. present(x1_max_slopes)) then
1407  sll_error("sll_s_cubic_spline_2d_init","hermite-periodic-case, it is not allowed to specify simultaneously a constant value for the slopes at x1_max and an array-specified set of slopes")
1408  end if
1409 
1410  ! But X1 slopes are.
1411  sll_allocate(this%x1_min_slopes(num_pts_x2), ierr)
1412  sll_allocate(this%x1_max_slopes(num_pts_x2), ierr)
1413 
1414  ! And since the X1 direction splines are computed second, then we
1415  ! need to convert the slope information into spline coefficient
1416  ! information.
1417  sll_allocate(this%x1_min_slopes_coeffs(0:num_pts_x2 + 2), ierr)
1418  sll_allocate(this%x1_max_slopes_coeffs(0:num_pts_x2 + 2), ierr)
1419 
1420  ! NOTE: because we are using the spline coefficients directly and not
1421  ! the slope values, the slope values are redundant and at this point
1422  ! we could deallocate those arrays...
1423 
1424  ! The following macro is obviously intended only for use within
1425  ! sll_s_cubic_spline_2d_init(). But this should be replaced with a subroutine
1426  ! whenever possible.
1427 
1428 #ifndef DOXYGEN_SHOULD_SKIP_THIS
1429 
1430 #define FILL_SLOPES(const_opt, input_opt, numpts, output, slopes) \
1431  if (present(input_opt)) then; \
1432  this%output(1:numpts) = input_opt(1:numpts); \
1433  this%slopes = .false.; \
1434  else if (present(const_opt)) then; \
1435  this%output(1:numpts) = const_opt; \
1436  this%slopes = .false.; \
1437  else; \
1438  this%slopes = .true.; \
1439  end if
1440 
1441 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
1442 
1443  ! Set the values of the slopes at x1_min
1444  fill_slopes(const_slope_x1_min, x1_min_slopes, num_pts_x2, x1_min_slopes, compute_slopes_x1_min)
1445 
1446  ! Set the values of the slopes at x1_max
1447  fill_slopes(const_slope_x1_max, x1_max_slopes, num_pts_x2, x1_max_slopes, compute_slopes_x1_max)
1448 
1449  case (9)
1450  ! Periodic in X1 and Hermite in X2
1451  if ( &
1452  present(x1_min_slopes) .or. present(x1_max_slopes) .or. &
1453  present(const_slope_x1_min) .or. present(const_slope_x1_max)) then
1454  sll_error("sll_s_cubic_spline_2d_init","periodic-hermite case, it is not allowed to specify the end slopes in the case of periodic boundary conditions.")
1455  sll_warning('sll_s_cubic_spline_2d_init', 'values of slopes in x1 are not taken into account as we are periodic in x1')
1456 
1457  end if
1458  if (present(const_slope_x2_min) .and. present(x2_min_slopes)) then
1459  sll_error("sll_s_cubic_spline_2d_init","periodic-hermite case, it is not allowed to specify simultaneously a constant value for the slopes at x2_min and an array-specified set of slopes")
1460  end if
1461  if (present(const_slope_x2_max) .and. present(x2_max_slopes)) then
1462  sll_error("sll_s_cubic_spline_2d_init","periodic-hermite case, it is not allowed to specify simultaneously a constant value for the slopes at x2_max and an array-specified set of slopes")
1463  end if
1464 
1465  ! But X2 slopes are.
1466  sll_allocate(this%x2_min_slopes(num_pts_x1), ierr)
1467  sll_allocate(this%x2_max_slopes(num_pts_x1), ierr)
1468  ! Except that the slope information is used directly, and not as
1469  ! spline coefficients data.
1470  ! Set the values of the slopes at x2_min
1471  fill_slopes(const_slope_x2_min, x2_min_slopes, num_pts_x1, x2_min_slopes, compute_slopes_x2_min)
1472 
1473  ! Set the values of the slopes at x2_max
1474  fill_slopes(const_slope_x2_max, x2_max_slopes, num_pts_x1, x2_max_slopes, compute_slopes_x2_max)
1475 
1476  case (10)
1477  ! Hermite conditions in both, X1 and X2
1478  if (present(const_slope_x1_min) .and. present(x1_min_slopes)) then
1479  sll_error("sll_s_cubic_spline_2d_init","hermite-hermite case, it is not allowed to specify simultaneously a constant value for the slopes at x1_min and an array-specified set of slopes")
1480  end if
1481  if (present(const_slope_x1_max) .and. present(x1_max_slopes)) then
1482  sll_error("sll_s_cubic_spline_2d_init","hermite-hermite case, it is not allowed to specify simultaneously a constant value for the slopes at x2_max and an array-specified set of slopes")
1483  end if
1484  if (present(const_slope_x2_min) .and. present(x2_min_slopes)) then
1485  sll_error("sll_s_cubic_spline_2d_init","hermite-hermite case, it is not allowed to specify simultaneously a constant value for the slopes at x2_min and an array-specified set of slopes")
1486  end if
1487  if (present(const_slope_x2_max) .and. present(x2_max_slopes)) then
1488  sll_error("sll_s_cubic_spline_2d_init","hermite-hermite case, it is not allowed to specify simultaneously a constant value for the slopes at x2_max and an array-specified set of slopes")
1489  end if
1490 
1491  ! Both, X1 and X2 slope arrays are needed.
1492  sll_allocate(this%x1_min_slopes(num_pts_x2), ierr)
1493  sll_allocate(this%x1_max_slopes(num_pts_x2), ierr)
1494  sll_allocate(this%x2_min_slopes(num_pts_x1), ierr)
1495  sll_allocate(this%x2_max_slopes(num_pts_x1), ierr)
1496  ! Given the order in which we compute the splines, first along X2 and
1497  ! then along X1, the coefficients-as-slopes data would never be used
1498  ! in the X2 direction. Consider eliminating these pointers from the
1499  ! object.
1500  sll_allocate(this%x1_min_slopes_coeffs(0:num_pts_x2 + 2), ierr)
1501  sll_allocate(this%x1_max_slopes_coeffs(0:num_pts_x2 + 2), ierr)
1502 
1503  ! Set the values of the slopes at x1_min
1504  fill_slopes(const_slope_x1_min, x1_min_slopes, num_pts_x2, x1_min_slopes, compute_slopes_x1_min)
1505 
1506  ! Set the values of the slopes at x1_max
1507  fill_slopes(const_slope_x1_max, x1_max_slopes, num_pts_x2, x1_max_slopes, compute_slopes_x1_max)
1508 
1509  ! Set the values of the slopes at x2_min
1510  fill_slopes(const_slope_x2_min, x2_min_slopes, num_pts_x1, x2_min_slopes, compute_slopes_x2_min)
1511 
1512  ! Set the values of the slopes at x2_max
1513  fill_slopes(const_slope_x2_max, x2_max_slopes, num_pts_x1, x2_max_slopes, compute_slopes_x2_max)
1514 
1515  case default
1516  sll_error("sll_s_cubic_spline_2d_init", "did not recognize given boundary conditions")
1517  end select
1518  ! Note: The indexing of the coefficients array includes the end-
1519  ! points 0, num_points, num_points+1, num_points+2. These are meant to
1520  ! store the boundary condition-specific data. The 'periodic' BC does
1521  ! not use the num_points+2 point.
1522  sll_allocate(this%coeffs(0:num_pts_x1 + 2, 0:num_pts_x2 + 2), ierr)
1523 
1524  end subroutine sll_s_cubic_spline_2d_init
1525 
1526  subroutine compute_spline_2d_prdc_prdc(data, spline)
1527  sll_real64, dimension(:, :), intent(in), target :: data ! data to be fit
1528  type(sll_t_cubic_spline_2d) :: spline
1529  sll_real64, dimension(:), pointer :: coeffs
1530  sll_int32 :: npx1
1531  sll_int32 :: npx2
1532  sll_real64, dimension(:), pointer :: d1
1533  sll_real64, dimension(:), pointer :: d2
1534  sll_real64, dimension(:), pointer :: datap ! 1D data slice pointer
1535  sll_int32 :: i
1536  sll_int32 :: j
1537  !if( .not. associated(spline) ) then
1538  ! ! FIXME: THROW ERROR
1539  ! print *, 'ERROR: compute_spline_2D_prdc_prdc(): ', &
1540  ! 'uninitialized spline object passed as argument. Exiting... '
1541  ! STOP
1542  !end if
1543  if ((size(data, 1) .lt. spline%num_pts_x1)) then
1544  ! FIXME: THROW ERROR
1545  print *, 'ERROR: compute_spline_2D_prdc_prdc(): '
1546  write (*, '(a, i8, a, i8)') 'spline object needs data of size >= ', &
1547  spline%num_pts_x1, ' . Passed size: ', size(data, 1)
1548  stop
1549  end if
1550  if ((size(data, 2) .lt. spline%num_pts_x2)) then
1551  ! FIXME: THROW ERROR
1552  print *, 'ERROR: compute_spline_2D_prdc_prdc(): '
1553  write (*, '(a, i8, a, i8)') 'spline object needs data of size >= ', &
1554  spline%num_pts_x2, ' . Passed size: ', size(data, 2)
1555  stop
1556  end if
1557  npx1 = spline%num_pts_x1
1558  npx2 = spline%num_pts_x2
1559  d1 => spline%d1
1560  d2 => spline%d2
1561  ! build splines along the x2 direction. Note: due to Fortran's
1562  ! column-major ordering, this uses long strides in memory.
1563  do i = 1, npx1
1564  datap => data(i, 1:npx2)
1565  coeffs => spline%coeffs(i, 0:npx2 + 2)
1566  call compute_spline_1d_periodic_aux(datap, npx2, spline%d2, coeffs)
1567  end do
1568  ! build splines along the x1 direction. Note: due to Fortran's
1569  ! column-major ordering, this involves short strides in memory.
1570  ! Note that the data are the spline coefficients computed in the
1571  ! previous step, so the array dimensions are slightly bigger than in
1572  ! the original data.
1573  do j = 0, npx2 + 2
1574  datap => spline%coeffs(1:npx1, j)
1575  ! same trick regarding the starting point of this pointer. This is
1576  ! not good.
1577  coeffs => spline%coeffs(0:npx1 + 2, j)
1578  call compute_spline_1d_periodic_aux(datap, npx1, d1, coeffs)
1579  end do
1580  end subroutine compute_spline_2d_prdc_prdc
1581 
1582  subroutine compute_spline_2d_hrmt_prdc(data, spline)
1583  sll_real64, dimension(:, :), intent(in), target :: data ! data to be fit
1584  type(sll_t_cubic_spline_2d) :: spline
1585  sll_real64, dimension(:), pointer :: coeffs
1586  sll_int32 :: npx1
1587  sll_int32 :: npx2
1588  sll_real64, dimension(:), pointer :: d1
1589  sll_real64, dimension(:), pointer :: d2
1590  sll_real64, dimension(:), pointer :: datap ! 1D data slice pointer
1591  sll_real64, dimension(:), pointer :: coeffs_ptr1
1592  sll_real64, dimension(:), pointer :: coeffs_ptr2
1593  sll_real64 :: min_slope ! slopes at endpoints
1594  sll_real64 :: max_slope
1595  sll_int32 :: i
1596  sll_int32 :: j
1597  sll_real64 :: r_x1_delta
1598 
1599  !if( .not. associated(spline) ) then
1600  ! ! FIXME: THROW ERROR
1601  ! print *, 'ERROR: compute_spline_2D_prdc_prdc(): ', &
1602  ! 'uninitialized spline object passed as argument. Exiting... '
1603  ! STOP
1604  !end if
1605  if ((size(data, 1) .lt. spline%num_pts_x1)) then
1606  ! FIXME: THROW ERROR
1607  print *, 'ERROR: compute_spline_2D_prdc_prdc(): '
1608  write (*, '(a, i8, a, i8)') 'spline object needs data of size >= ', &
1609  spline%num_pts_x1, ' . Passed size: ', size(data, 1)
1610  stop
1611  end if
1612  if ((size(data, 2) .lt. spline%num_pts_x2)) then
1613  ! FIXME: THROW ERROR
1614  print *, 'ERROR: compute_spline_2D_prdc_prdc(): '
1615  write (*, '(a, i8, a, i8)') 'spline object needs data of size >= ', &
1616  spline%num_pts_x2, ' . Passed size: ', size(data, 2)
1617  stop
1618  end if
1619  npx1 = spline%num_pts_x1
1620  npx2 = spline%num_pts_x2
1621  d1 => spline%d1
1622  d2 => spline%d2
1623  ! build splines along the x2 direction (the periodic direction).
1624  ! Note: due to Fortran's column-major ordering, this uses long strides
1625  ! in memory.
1626  do i = 1, npx1
1627  datap => data(i, 1:npx2)
1628  coeffs => spline%coeffs(i, 0:npx2 + 2)
1629  call compute_spline_1d_periodic_aux(datap, npx2, d2, coeffs)
1630  end do
1631  ! build splines along the x1 direction (the Hermite direction).
1632  ! Note: due to Fortran's column-major ordering, this involves short
1633  ! strides in memory.
1634  !
1635  ! Note also that the data are the spline coefficients computed in the
1636  ! previous step, so the array dimensions are slightly bigger than in
1637  ! the original data.
1638 
1639  ! First, compute the spline coefficients along the lines (:,0) and
1640  ! (:,npx2+1)
1641  r_x1_delta = 1.0_f64/spline%x1_delta
1642  if (spline%compute_slopes_x1_min .eqv. .true.) then
1643  ! compute default value for the derivative at the first point based
1644  ! on the given data. Can't use the macros here... go figure...
1645  do j = 1, npx2
1646  spline%x1_min_slopes(j) = &
1647  r_x1_delta*(-(25.0_f64/12.0_f64)*data(1, j) + &
1648  4.0_f64*data(2, j) - &
1649  3.0_f64*data(3, j) + &
1650  (4.0_f64/3.0_f64)*data(4, j) - &
1651  0.25_f64*data(5, j))
1652  end do
1653  end if
1654  if (spline%compute_slopes_x1_max .eqv. .true.) then
1655  ! estimate the derivative at the last point.
1656  do j = 1, npx2
1657  spline%x1_max_slopes(j) = &
1658  r_x1_delta*(+0.25_f64*data(npx1 - 4, j) - &
1659  (4.0_f64/3.0_f64)*data(npx1 - 3, j) + &
1660  3.0_f64*data(npx1 - 2, j) - &
1661  4.0_f64*data(npx1 - 1, j) + &
1662  (25.0_f64/12.0_f64)*data(npx1, j))
1663  end do
1664  end if
1665  ! At this point, the values of the slopes are available because either
1666  ! the caller has supplied the values or they have been computed
1667  ! numerically.
1668  ! Compute the spline coefficients for the available slope values
1669  coeffs_ptr1 => spline%x1_min_slopes_coeffs(0:npx2 + 2)
1671  spline%x1_min_slopes, &
1672  npx2, &
1673  spline%d2, &
1674  coeffs_ptr1)
1675  coeffs_ptr2 => spline%x1_max_slopes_coeffs(0:npx2 + 2)
1677  spline%x1_max_slopes, &
1678  npx2, &
1679  spline%d2, &
1680  coeffs_ptr2)
1681 
1682  do j = 0, npx2 + 2
1683  datap => spline%coeffs(1:npx1, j)
1684  coeffs => spline%coeffs(0:npx1 + 2, j)
1685  min_slope = spline%x1_min_slopes_coeffs(j)
1686  max_slope = spline%x1_max_slopes_coeffs(j)
1688  datap, &
1689  npx1, &
1690  spline%d1, &
1691  min_slope, &
1692  max_slope, &
1693  spline%x1_delta, &
1694  coeffs)
1695 ! call compute_spline_1D_periodic_aux( datap, npx1, d1, coeffs )
1696  end do
1697  end subroutine compute_spline_2d_hrmt_prdc
1698 
1699  subroutine compute_spline_2d_prdc_hrmt(data, spline)
1700  sll_real64, dimension(:, :), intent(in), target :: data ! data to be fit
1701  type(sll_t_cubic_spline_2d) :: spline
1702  sll_real64, dimension(:), pointer :: coeffs
1703  sll_int32 :: npx1
1704  sll_int32 :: npx2
1705  sll_real64, dimension(:), pointer :: d1
1706  sll_real64, dimension(:), pointer :: d2
1707  sll_real64, dimension(:), pointer :: datap ! 1D data slice pointer
1708  sll_real64 :: min_slope ! slopes at endpoints
1709  sll_real64 :: max_slope
1710  sll_int32 :: i
1711  sll_int32 :: j
1712  sll_real64 :: r_x2_delta
1713 
1714  !if( .not. associated(spline) ) then
1715  ! ! FIXME: THROW ERROR
1716  ! print *, 'ERROR: compute_spline_2D_prdc_prdc(): ', &
1717  ! 'uninitialized spline object passed as argument. Exiting... '
1718  ! STOP
1719  !end if
1720  if ((size(data, 1) .lt. spline%num_pts_x1)) then
1721  ! FIXME: THROW ERROR
1722  print *, 'ERROR: compute_spline_2D_prdc_prdc(): '
1723  write (*, '(a, i8, a, i8)') 'spline object needs data of size >= ', &
1724  spline%num_pts_x1, ' . Passed size: ', size(data, 1)
1725  stop
1726  end if
1727  if ((size(data, 2) .lt. spline%num_pts_x2)) then
1728  ! FIXME: THROW ERROR
1729  print *, 'ERROR: compute_spline_2D_prdc_prdc(): '
1730  write (*, '(a, i8, a, i8)') 'spline object needs data of size >= ', &
1731  spline%num_pts_x2, ' . Passed size: ', size(data, 2)
1732  stop
1733  end if
1734  npx1 = spline%num_pts_x1
1735  npx2 = spline%num_pts_x2
1736  d1 => spline%d1
1737  d2 => spline%d2
1738  r_x2_delta = 1.0_f64/spline%x2_delta
1739 
1740  ! Numerically compute the values of the slopes if the user has not given
1741  ! them.
1742  if (spline%compute_slopes_x2_min .eqv. .true.) then
1743  ! forward difference scheme to estimate the derivative
1744  ! at the first point based on the given data.
1745  do i = 1, npx1
1746  spline%x2_min_slopes(i) = &
1747  r_x2_delta*(-(25.0_f64/12.0_f64)*data(i, 1) + &
1748  4.0_f64*data(i, 2) - &
1749  3.0_f64*data(i, 3) + &
1750  (4.0_f64/3.0_f64)*data(i, 4) - &
1751  0.25_f64*data(i, 5))
1752  end do
1753  end if
1754  if (spline%compute_slopes_x2_max .eqv. .true.) then
1755  ! backward difference scheme to estimate the derivative
1756  ! at the last point.
1757  do i = 1, npx1
1758  spline%x2_max_slopes(i) = &
1759  r_x2_delta*(+0.25_f64*data(i, npx2 - 4) - &
1760  (4.0_f64/3.0_f64)*data(i, npx2 - 3) + &
1761  3.0_f64*data(i, npx2 - 2) - &
1762  4.0_f64*data(i, npx2 - 1) + &
1763  (25.0_f64/12.0_f64)*data(i, npx2))
1764  end do
1765  end if
1766 
1767  ! build splines along the x2 direction (hermite direction). Note: due
1768  ! to Fortran's column-major ordering, this uses long strides in memory.
1769  do i = 1, npx1
1770  datap => data(i, 1:npx2)
1771  coeffs => spline%coeffs(i, 0:npx2 + 2)
1772  min_slope = spline%x2_min_slopes(i)
1773  max_slope = spline%x2_max_slopes(i)
1775  datap, &
1776  npx2, &
1777  spline%d2, &
1778  min_slope, &
1779  max_slope, &
1780  spline%x2_delta, &
1781  coeffs)
1782  end do
1783  ! build splines along the x1 direction. Note: due to Fortran's
1784  ! column-major ordering, this involves short strides in memory.
1785  ! Note that the data are the spline coefficients computed in the
1786  ! previous step, so the array dimensions are slightly bigger than in
1787  ! the original data.
1788  do j = 0, npx2 + 1
1789  datap => spline%coeffs(1:npx1, j)
1790  coeffs => spline%coeffs(0:npx1 + 2, j)
1791  call compute_spline_1d_periodic_aux(datap, npx1, spline%d1, coeffs)
1792  end do
1793  end subroutine compute_spline_2d_prdc_hrmt
1794 
1795  subroutine compute_spline_2d_hrmt_hrmt(data, spline)
1796  sll_real64, dimension(:, :), intent(in), target :: data ! data to be fit
1797  type(sll_t_cubic_spline_2d) :: spline
1798  sll_real64, dimension(:), pointer :: coeffs
1799  sll_int32 :: npx1
1800  sll_int32 :: npx2
1801  sll_real64, dimension(:), pointer :: d1
1802  sll_real64, dimension(:), pointer :: d2
1803  sll_real64, dimension(:), pointer :: datap ! 1D data slice pointer
1804  sll_real64 :: min_slope ! slopes at endpoints
1805  sll_real64 :: max_slope
1806  sll_int32 :: i
1807  sll_int32 :: j
1808  sll_real64 :: r_x1_delta ! reciprocal of x1_delta
1809  sll_real64 :: r_x2_delta ! reciprocal of x2_delta
1810 
1811  !if( .not. associated(spline) ) then
1812  ! ! FIXME: THROW ERROR
1813  ! print *, 'ERROR: compute_spline_2D_hrmt_hrmt(): ', &
1814  ! 'uninitialized spline object passed as argument. Exiting... '
1815  ! STOP
1816  !end if
1817  if ((size(data, 1) .lt. spline%num_pts_x1)) then
1818  ! FIXME: THROW ERROR
1819  print *, 'ERROR: compute_spline_2D_hrmt_hrmt(): '
1820  write (*, '(a, i8, a, i8)') 'spline object needs data of size >= ', &
1821  spline%num_pts_x1, ' . Passed size: ', size(data, 1)
1822  stop
1823  end if
1824  if ((size(data, 2) .lt. spline%num_pts_x2)) then
1825  ! FIXME: THROW ERROR
1826  print *, 'ERROR: compute_spline_2D_hrmt_hrmt(): '
1827  write (*, '(a, i8, a, i8)') 'spline object needs data of size >= ', &
1828  spline%num_pts_x2, ' . Passed size: ', size(data, 2)
1829  stop
1830  end if
1831  npx1 = spline%num_pts_x1
1832  npx2 = spline%num_pts_x2
1833  d1 => spline%d1
1834  d2 => spline%d2
1835  r_x1_delta = 1.0_f64/spline%x1_delta
1836  r_x2_delta = 1.0_f64/spline%x2_delta
1837 
1838  ! Compute the values of the slopes in case that they were not given.
1839  if (spline%compute_slopes_x1_min .eqv. .true.) then
1840  ! forward difference scheme to estimate the derivative
1841  ! at the first point based on the given data.
1842  do j = 1, npx2
1843  spline%x1_min_slopes(j) = &
1844  r_x1_delta*(-(25.0_f64/12.0_f64)*data(1, j) + &
1845  4.0_f64*data(2, j) - &
1846  3.0_f64*data(3, j) + &
1847  (4.0_f64/3.0_f64)*data(4, j) - &
1848  0.25_f64*data(5, j))
1849  end do
1850  end if
1851  if (spline%compute_slopes_x1_max .eqv. .true.) then
1852  ! backward difference scheme to estimate the derivative
1853  ! at the last point.
1854  do j = 1, npx2
1855  spline%x1_max_slopes(j) = &
1856  r_x1_delta*(+0.25_f64*data(npx1 - 4, j) - &
1857  (4.0_f64/3.0_f64)*data(npx1 - 3, j) + &
1858  3.0_f64*data(npx1 - 2, j) - &
1859  4.0_f64*data(npx1 - 1, j) + &
1860  (25.0_f64/12.0_f64)*data(npx1, j))
1861  end do
1862  end if
1863 
1864  if (spline%compute_slopes_x2_min .eqv. .true.) then
1865  ! forward difference scheme to estimate the derivative
1866  ! at the first point based on the given data.
1867  do i = 1, npx1
1868  spline%x2_min_slopes(i) = &
1869  r_x2_delta*(-(25.0_f64/12.0_f64)*data(i, 1) + &
1870  4.0_f64*data(i, 2) - &
1871  3.0_f64*data(i, 3) + &
1872  (4.0_f64/3.0_f64)*data(i, 4) - &
1873  0.25_f64*data(i, 5))
1874  end do
1875  end if
1876  if (spline%compute_slopes_x2_max .eqv. .true.) then
1877  ! backward difference scheme to estimate the derivative
1878  ! at the last point.
1879  do i = 1, npx1
1880  spline%x2_max_slopes(i) = &
1881  r_x2_delta*(+0.25_f64*data(i, npx2 - 4) - &
1882  (4.0_f64/3.0_f64)*data(i, npx2 - 3) + &
1883  3.0_f64*data(i, npx2 - 2) - &
1884  4.0_f64*data(i, npx2 - 1) + &
1885  (25.0_f64/12.0_f64)*data(i, npx2))
1886  end do
1887  end if
1888 
1889  ! build splines along the x2 direction. Note: due to Fortran's
1890  ! column-major ordering, this uses long strides in memory.
1891  do i = 1, npx1
1892  datap => data(i, 1:npx2)
1893  coeffs => spline%coeffs(i, 0:npx2 + 2)
1894  min_slope = spline%x2_min_slopes(i)
1895  max_slope = spline%x2_max_slopes(i)
1897  datap, &
1898  npx2, &
1899  spline%d2, &
1900  min_slope, &
1901  max_slope, &
1902  spline%x2_delta, &
1903  coeffs)
1904  end do
1905  ! build splines along the x1 direction. Note: due to Fortran's
1906  ! column-major ordering, this involves short strides in memory.
1907  ! Note that the data are the spline coefficients computed in the
1908  ! previous step, so the array dimensions are slightly bigger than in
1909  ! the original data.
1910 
1911  ! First we compute out of the loop the splines for the coefficients
1912  ! in the (:,0) and (:,npx2+1) rows. TO DO THIS PROPERLY WE SHOULD
1913  ! INTRODUCE AN ESTIMATE OF THE SLOPES, WHICH FOR THESE VALUES FALL
1914  ! OUT OF RANGE. Here we just "reflect" the values around the edge point.
1915  ! ... it's better than nothing.
1916  datap => spline%coeffs(1:npx1, 0)
1917  coeffs => spline%coeffs(0:npx1 + 2, 0)
1918  min_slope = spline%x1_min_slopes(2)
1919  max_slope = spline%x1_max_slopes(2)
1921  datap, &
1922  npx1, &
1923  spline%d1, &
1924  min_slope, &
1925  max_slope, &
1926  spline%x1_delta, &
1927  coeffs)
1928 
1929  datap => spline%coeffs(1:npx1, npx2 + 1)
1930  coeffs => spline%coeffs(0:npx1 + 2, npx2 + 1)
1931  min_slope = spline%x1_min_slopes(npx2 - 1)
1932  max_slope = spline%x1_max_slopes(npx2 - 1)
1934  datap, &
1935  npx1, &
1936  spline%d1, &
1937  min_slope, &
1938  max_slope, &
1939  spline%x1_delta, &
1940  coeffs)
1941 
1942  do j = 1, npx2
1943  datap => spline%coeffs(1:npx1, j)
1944  coeffs => spline%coeffs(0:npx1 + 2, j)
1945  min_slope = spline%x1_min_slopes(j)
1946  max_slope = spline%x1_max_slopes(j)
1948  datap, &
1949  npx1, &
1950  spline%d1, &
1951  min_slope, &
1952  max_slope, &
1953  spline%x1_delta, &
1954  coeffs)
1955  end do
1956  end subroutine compute_spline_2d_hrmt_hrmt
1957 
1965  sll_real64, dimension(:, :), intent(in), target :: data ! data to be fit
1966  type(sll_t_cubic_spline_2d) :: spline
1967  sll_int32 :: bc1
1968  sll_int32 :: bc2
1969  sll_int32 :: bc_selector
1970 
1971  !if( .not. associated(spline) ) then
1972  ! ! FIXME: THROW ERROR
1973  ! print *, 'ERROR: compute_spline_2D(): ', &
1974  ! 'uninitialized spline object passed as argument. Exiting... '
1975  ! STOP
1976  !end if
1977 
1978  bc1 = spline%x1_bc_type
1979  bc2 = spline%x2_bc_type
1980 
1981  ! Treat the bc_selector variable essentially like a bit field, to
1982  ! accumulate the information on the different boundary conditions
1983  ! given. This scheme allows to add more types of boundary conditions
1984  ! if necessary.
1985  bc_selector = 0
1986 
1987  ! We make every case explicit to facilitate adding more BC types in
1988  ! the future.
1989  if (bc1 .eq. sll_p_periodic) then
1990  bc_selector = bc_selector + 1
1991  end if
1992 
1993  if (bc1 .eq. sll_p_hermite) then
1994  bc_selector = bc_selector + 2
1995  end if
1996 
1997  if (bc2 .eq. sll_p_periodic) then
1998  bc_selector = bc_selector + 4
1999  end if
2000 
2001  if (bc2 .eq. sll_p_hermite) then
2002  bc_selector = bc_selector + 8
2003  end if
2004 
2005  select case (bc_selector)
2006  case (5)
2007  ! both boundary condition types are periodic
2008  call compute_spline_2d_prdc_prdc(data, spline)
2009  case (6)
2010  ! hermite in X1 and periodic in X2
2011  call compute_spline_2d_hrmt_prdc(data, spline)
2012  case (9)
2013  ! periodic condition in X1 and hermite in X2
2014  call compute_spline_2d_prdc_hrmt(data, spline)
2015  case (10)
2016  ! Hermite conditions in both, X1 and X2
2017  call compute_spline_2d_hrmt_hrmt(data, spline)
2018  case default
2019  print *, 'ERROR: compute_spline_2D(): ', &
2020  'did not recognize given boundary condition combination.'
2021  stop
2022  end select
2024 
2036  subroutine sll_s_cubic_spline_2d_deposit_value(x1, x2, spline, a_out)
2037  intrinsic :: real, int
2038  sll_real64, dimension(1:, 1:), intent(in) :: x1
2039  sll_real64, dimension(1:, 1:), intent(in) :: x2
2040  type(sll_t_cubic_spline_2d) :: spline
2041  sll_real64, dimension(:, :), intent(out) :: a_out
2042 
2043  sll_real64 :: cij ! C_ij
2044  sll_real64 :: x1_min
2045  sll_real64 :: x2_min
2046  sll_real64 :: dx1
2047  sll_real64 :: dx2
2048  sll_real64 :: cdx1 ! 1-dx1
2049  sll_real64 :: cdx2 ! 1-dx2
2050  sll_int32 :: cell1
2051  sll_int32 :: cell2
2052  sll_real64 :: rh1
2053  sll_real64 :: rh2
2054  sll_int32 :: n1
2055  sll_int32 :: n2
2056 
2057  ! local variables
2058  sll_int32 :: i1
2059  sll_int32 :: i2
2060 
2061  sll_int32 :: nt1
2062  sll_int32 :: nt2
2063 
2064  sll_real64 :: svalx1, svalx2, svalx3, svalx4
2065  sll_real64 :: svaly1, svaly2, svaly3, svaly4
2066  sll_real64 :: ax1, ax2, ax3, ay1, ay2, ay3
2067 
2068  sll_real64 :: t1, t2
2069  sll_int32 :: ipm1, ip, ipp1, ipp2
2070  sll_int32 :: jpm1, jp, jpp1, jpp2
2071 
2072  sll_int32 :: bc1
2073  sll_int32 :: bc2
2074 
2075  !if( .not. associated(spline) ) then
2076  ! ! FIXME: THROW ERROR
2077  ! print *, 'ERROR: sll_s_cubic_spline_2d_deposit_value(): ', &
2078  ! 'uninitialized spline object passed as argument. Exiting... '
2079  ! STOP
2080  !end if
2081 
2082  if ((size(x1, 1) .ne. spline%num_pts_x1) .or. (size(x1, 2) .ne. spline%num_pts_x2)) then
2083  ! FIXME: THROW ERROR
2084  print *, 'ERROR: sll_s_cubic_spline_2d_deposit_value(): '
2085  write (*, '(a, i8, i8, a, i8, i8, a)') 'array of feets of characteristics needs data of size = (', &
2086  spline%num_pts_x1, spline%num_pts_x2, ') . Passed size: (', size(x1, 1), size(x1, 2), ')'
2087  stop
2088  end if
2089  if ((size(x2, 1) .ne. spline%num_pts_x1) .or. (size(x2, 2) .ne. spline%num_pts_x2)) then
2090  ! FIXME: THROW ERROR
2091  print *, 'ERROR: sll_s_cubic_spline_2d_deposit_value(): '
2092  write (*, '(a, i8, i8, a, i8, i8, a)') 'array of feets of characteristics needs data of size = (', &
2093  spline%num_pts_x1, spline%num_pts_x2, ') . Passed size: (', size(x2, 1), size(x2, 2), ')'
2094  stop
2095  end if
2096  if ((size(a_out, 1) .ne. spline%num_pts_x1) .or. (size(a_out, 2) .ne. spline%num_pts_x2)) then
2097  ! FIXME: THROW ERROR
2098  print *, 'ERROR: sll_s_cubic_spline_2d_deposit_value(): '
2099  write (*, '(a, i8, i8, a, i8, i8, a)') 'array of feets of characteristics needs data of size = (', &
2100  spline%num_pts_x1, spline%num_pts_x2, ') . Passed size: (', size(a_out, 1), size(a_out, 2), ')'
2101  stop
2102  end if
2103 
2104  bc1 = spline%x1_bc_type
2105  bc2 = spline%x2_bc_type
2106 
2107  x1_min = spline%x1_min
2108  x2_min = spline%x2_min
2109  rh1 = spline%x1_rdelta
2110  rh2 = spline%x2_rdelta
2111 
2112  n1 = spline%num_pts_x1
2113  n2 = spline%num_pts_x2
2114 
2115  if (bc1 .eq. sll_p_periodic) then
2116  nt1 = n1 - 1
2117  end if
2118  if (bc1 .eq. sll_p_hermite) then
2119  nt1 = n1
2120  end if
2121  if (bc2 .eq. sll_p_periodic) then
2122  nt2 = n2 - 1
2123  end if
2124  if (bc2 .eq. sll_p_hermite) then
2125  nt2 = n2
2126  end if
2127 
2128  a_out = 0._f64
2129 
2130  do i1 = 1, nt1
2131  do i2 = 1, nt2
2132 
2133  ! find the cell and offset for x1
2134  t1 = (x1(i1, i2) - x1_min)*rh1
2135  cell1 = floor(t1) + 1
2136  dx1 = t1 - real(cell1 - 1, f64)
2137  cdx1 = 1.0_f64 - dx1
2138 
2139  ! find the cell and offset for x2
2140  t2 = (x2(i1, i2) - x2_min)*rh2
2141  cell2 = floor(t2) + 1
2142  dx2 = t2 - real(cell2 - 1, f64)
2143  cdx2 = 1.0_f64 - dx2
2144 
2145  ! checking if the particle is well localized
2146  ! the particle is in the cell (cell1,cell2)
2147  ! placed in ((cell1-1)/rh1+x1_min,(cell2-1)/rh2+x2_min)
2148  if ((real(cell1 - 1, f64)/rh1 + x1_min > x1(i1, i2)) .or. (x1(i1, i2) > real(cell1, f64)/rh1 + x1_min)) then
2149  print *, 'problem with the localization of r', &
2150  real(cell1 - 1, f64)/rh1 + x1_min, x1(i1, i2), real(cell1, f64)/rh1 + x1_min, dx1
2151  end if
2152  if ((real(cell2 - 1, f64)/rh2 + x2_min > x2(i1, i2)) .or. (x2(i1, i2) > real(cell2, f64)/rh2 + x2_min)) then
2153  print *, 'problem with the localization of theta', &
2154  real(cell2 - 1, f64)/rh2 + x2_min, x2(i1, i2), real(cell2, f64)/rh2 + x2_min, dx2
2155  end if
2156 
2157  cij = spline%coeffs(i1, i2)
2158 
2159  ! index depending on the BC
2160  if (bc1 .eq. sll_p_periodic) then
2161  ipm1 = mod(cell1 + n1 - 3, n1 - 1) + 1
2162  ip = mod(cell1 + n1 - 2, n1 - 1) + 1
2163  ipp1 = mod(cell1 + n1 - 1, n1 - 1) + 1
2164  ipp2 = mod(cell1 + n1, n1 - 1) + 1
2165  end if
2166  if (bc1 .eq. sll_p_hermite) then
2167  ipm1 = cell1 - 1
2168  ip = cell1
2169  ipp1 = cell1 + 1
2170  ipp2 = cell1 + 2
2171  end if
2172  if (bc2 .eq. sll_p_periodic) then
2173  jpm1 = mod(cell2 + n2 - 3, n2 - 1) + 1
2174  jp = mod(cell2 + n2 - 2, n2 - 1) + 1
2175  jpp1 = mod(cell2 + n2 - 1, n2 - 1) + 1
2176  jpp2 = mod(cell2 + n2, n2 - 1) + 1
2177  end if
2178  if (bc2 .eq. sll_p_hermite) then
2179  jpm1 = cell2 - 1
2180  jp = cell2
2181  jpp1 = cell2 + 1
2182  jpp2 = cell2 + 2
2183  end if
2184 
2185  ax1 = cdx1
2186  ax2 = ax1*ax1
2187  ax3 = ax2*ax1
2188 
2189  ay1 = cdx2
2190  ay2 = ay1*ay1
2191  ay3 = ay2*ay1
2192 
2193  svalx1 = ax3/6._f64
2194  svalx2 = (0.5_f64*(-ax3 + ax2 + ax1) + 1._f64/6._f64)
2195  svalx3 = (2._f64/3._f64 - ax2 + 0.5_f64*ax3)
2196  svalx4 = ((1._f64 - ax3)/6._f64 + 0.5_f64*(ax2 - ax1))
2197 
2198  svaly1 = ay3/6._f64
2199  svaly2 = 0.5_f64*(-ay3 + ay2 + ay1) + 1._f64/6._f64
2200  svaly3 = 2._f64/3.0_f64 - ay2 + 0.5_f64*ay3
2201  svaly4 = (1._f64 - ay3)/6._f64 + 0.5_f64*(ay2 - ay1)
2202 
2203  if (ipm1 .ge. 1) then
2204  if (jpm1 .ge. 1) then
2205  a_out(ipm1, jpm1) = a_out(ipm1, jpm1) + cij*svalx1*svaly1
2206  end if
2207  a_out(ipm1, jp) = a_out(ipm1, jp) + cij*svalx1*svaly2
2208  if (jpp1 .le. n2) then
2209  a_out(ipm1, jpp1) = a_out(ipm1, jpp1) + cij*svalx1*svaly3
2210  end if
2211  if (jpp2 .le. n2) then
2212  a_out(ipm1, jpp2) = a_out(ipm1, jpp2) + cij*svalx1*svaly4
2213  end if
2214  end if
2215 
2216  if (jpm1 .ge. 1) then
2217  a_out(ip, jpm1) = a_out(ip, jpm1) + cij*svalx2*svaly1
2218  end if
2219  a_out(ip, jp) = a_out(ip, jp) + cij*svalx2*svaly2
2220  if (jpp1 .le. n2) then
2221  a_out(ip, jpp1) = a_out(ip, jpp1) + cij*svalx2*svaly3
2222  end if
2223  if (jpp2 .le. n2) then
2224  a_out(ip, jpp2) = a_out(ip, jpp2) + cij*svalx2*svaly4
2225  end if
2226 
2227  if (ipp1 .le. n1) then
2228  if (jpm1 .ge. 1) then
2229  a_out(ipp1, jpm1) = a_out(ipp1, jpm1) + cij*svalx3*svaly1
2230  end if
2231  a_out(ipp1, jp) = a_out(ipp1, jp) + cij*svalx3*svaly2
2232  if (jpp1 .le. n2) then
2233  a_out(ipp1, jpp1) = a_out(ipp1, jpp1) + cij*svalx3*svaly3
2234  end if
2235  if (jpp2 .le. n2) then
2236  a_out(ipp1, jpp2) = a_out(ipp1, jpp2) + cij*svalx3*svaly4
2237  end if
2238  end if
2239 
2240  if (ipp2 .le. n1) then
2241  if (jpm1 .ge. 1) then
2242  a_out(ipp2, jpm1) = a_out(ipp2, jpm1) + cij*svalx4*svaly1
2243  end if
2244  a_out(ipp2, jp) = a_out(ipp2, jp) + cij*svalx4*svaly2
2245  if (jpp1 .le. n2) then
2246  a_out(ipp2, jpp1) = a_out(ipp2, jpp1) + cij*svalx4*svaly3
2247  end if
2248  if (jpp2 .le. n2) then
2249  a_out(ipp2, jpp2) = a_out(ipp2, jpp2) + cij*svalx4*svaly4
2250  end if
2251  end if
2252 
2253  if (bc1 .eq. sll_p_hermite) then
2254  if (i1 .eq. 1) then
2255  if (jpm1 .ge. 1) then
2256  a_out(1, jpm1) = a_out(1, jpm1) + spline%coeffs(0, i2)/6._f64*svaly1
2257  end if
2258  a_out(1, jp) = a_out(1, jp) + spline%coeffs(0, i2)/6._f64*svaly2
2259  if (jpp1 .le. n2) then
2260  a_out(1, jpp1) = a_out(1, jpp1) + spline%coeffs(0, i2)/6._f64*svaly3
2261  end if
2262  if (jpp2 .le. n2) then
2263  a_out(1, jpp2) = a_out(1, jpp2) + spline%coeffs(0, i2)/6._f64*svaly4
2264  end if
2265  end if
2266 
2267  if (i1 .eq. n1) then
2268  if (jpm1 .ge. 1) then
2269  a_out(n1, jpm1) = a_out(n1, jpm1) + spline%coeffs(n1 + 1, i2)/6._f64*svaly1
2270  end if
2271  a_out(n1, jp) = a_out(n1, jp) + spline%coeffs(n1 + 1, i2)/6._f64*svaly2
2272  if (jpp1 .le. n2) then
2273  a_out(n1, jpp1) = a_out(n1, jpp1) + spline%coeffs(n1 + 1, i2)/6._f64*svaly3
2274  end if
2275  if (jpp2 .le. n2) then
2276  a_out(n1, jpp2) = a_out(n1, jpp2) + spline%coeffs(n1 + 1, i2)/6._f64*svaly4
2277  end if
2278  end if
2279  end if
2280 
2281  if (bc2 .eq. sll_p_hermite) then
2282  if (i2 .eq. 1) then
2283  if (ipm1 .ge. 1) then
2284  a_out(ipm1, 1) = a_out(ipm1, 1) + spline%coeffs(i1, 0)/6._f64*svalx1
2285  end if
2286  a_out(ip, 1) = a_out(ip, 1) + spline%coeffs(i1, 0)/6._f64*svalx2
2287  if (ipp1 .ne. n1) then
2288  a_out(ipp1, 1) = a_out(ipp1, 1) + spline%coeffs(i1, 0)/6._f64*svalx3
2289  end if
2290  if (ipp2 .ne. n1) then
2291  a_out(ipp2, 1) = a_out(ipp2, 1) + spline%coeffs(i1, 0)/6._f64*svalx4
2292  end if
2293  end if
2294 
2295  if (i2 .eq. n2) then
2296  if (ipm1 .ge. 1) then
2297  a_out(ipm1, n2) = a_out(ipm1, n2) + spline%coeffs(i1, n2 + 1)/6._f64*svalx1
2298  end if
2299  a_out(ip, n1) = a_out(ip, n1) + spline%coeffs(i1, n2 + 1)/6._f64*svalx2
2300  if (ipp1 .ne. n1) then
2301  a_out(ipp1, n2) = a_out(ipp1, n2) + spline%coeffs(i1, n2 + 1)/6._f64*svalx3
2302  end if
2303  if (ipp2 .ne. n1) then
2304  a_out(ipp2, n2) = a_out(ipp2, n2) + spline%coeffs(i1, n2 + 1)/6._f64*svalx4
2305  end if
2306  end if
2307  end if
2308 
2309  end do
2310  end do
2311 
2312  if (bc1 .eq. sll_p_periodic) then
2313  a_out(n1, :) = a_out(1, :)
2314  end if
2315  if (bc2 .eq. sll_p_periodic) then
2316  a_out(:, n2) = a_out(:, 1)
2317  end if
2318 
2320 
2328  function sll_f_cubic_spline_2d_eval(x1, x2, spline)
2329  sll_real64 :: sll_f_cubic_spline_2d_eval
2330  intrinsic :: associated, int, real
2331  sll_real64, intent(in) :: x1
2332  sll_real64, intent(in) :: x2
2333  type(sll_t_cubic_spline_2d), intent(in) :: spline
2334  sll_real64 :: rh1 ! reciprocal of cell spacing
2335  sll_real64 :: rh2 ! reciprocal of cell spacing
2336  sll_int32 :: cell
2337  sll_real64 :: dx
2338  sll_real64 :: cdx ! 1-dx
2339  sll_real64 :: t0 ! temp/scratch variables ...
2340  sll_real64 :: t1
2341  sll_real64 :: t2
2342  sll_real64 :: t3
2343  sll_real64 :: t4
2344  sll_real64 :: cim1 ! C_(i-1)
2345  sll_real64 :: ci ! C_i
2346  sll_real64 :: cip1 ! C_(i+1)
2347  sll_real64 :: cip2 ! C_(i+2)
2348  sll_real64 :: x1_min
2349  sll_real64 :: x2_min
2350  sll_int32 :: num_pts_x1
2351  sll_int32 :: num_pts_x2
2352  sll_real64, dimension(:), pointer :: coeffs_line_jm1
2353  sll_real64, dimension(:), pointer :: coeffs_line_j
2354  sll_real64, dimension(:), pointer :: coeffs_line_jp1
2355  sll_real64, dimension(:), pointer :: coeffs_line_jp2
2356  ! We set these as assertions since we want the flexibility of turning
2357  ! them off.
2358  sll_assert((x1 .ge. spline%x1_min) .and. (x1 .le. spline%x1_max))
2359  sll_assert((x2 .ge. spline%x2_min) .and. (x2 .le. spline%x2_max))
2360 
2361  x1_min = spline%x1_min
2362  x2_min = spline%x2_min
2363  num_pts_x1 = spline%num_pts_x1
2364  num_pts_x2 = spline%num_pts_x2
2365  rh1 = spline%x1_rdelta
2366  rh2 = spline%x2_rdelta
2367  ! find the cell and offset for x2
2368  t0 = (x2 - x2_min)*rh2
2369  cell = int(t0) + 1
2370  dx = t0 - real(cell - 1, f64)
2371  cdx = 1.0_f64 - dx
2372  ! write (*,'(a,i8, a, f20.12)') 'cell = ', cell, ', dx = ', dx
2373  ! interpolate the coefficients along the line of constant x1. These
2374  ! computations are independent from one another. A little problem is
2375  ! the redundancy in the computation of the cell and offset along each
2376  ! of the constant x2 lines, as this will be done by each call of
2377  ! interpolate_value_aux(). This suggests that the proper refactoring
2378  ! of this function would have the cell and offset as arguments.
2379  coeffs_line_jm1 => spline%coeffs(0:num_pts_x1 + 2, cell - 1)
2380  coeffs_line_j => spline%coeffs(0:num_pts_x1 + 2, cell)
2381  coeffs_line_jp1 => spline%coeffs(0:num_pts_x1 + 2, cell + 1)
2382  coeffs_line_jp2 => spline%coeffs(0:num_pts_x1 + 2, cell + 2)
2383  cim1 = interpolate_value_aux(x1, x1_min, rh1, coeffs_line_jm1)
2384  ci = interpolate_value_aux(x1, x1_min, rh1, coeffs_line_j)
2385  cip1 = interpolate_value_aux(x1, x1_min, rh1, coeffs_line_jp1)
2386  cip2 = interpolate_value_aux(x1, x1_min, rh1, coeffs_line_jp2)
2387  t1 = 3.0_f64*ci
2388  t3 = 3.0_f64*cip1
2389  t2 = cdx*(cdx*(cdx*(cim1 - t1) + t1) + t1) + ci
2390  t4 = dx*(dx*(dx*(cip2 - t3) + t3) + t3) + cip1
2391  sll_f_cubic_spline_2d_eval = (1.0_f64/6.0_f64)*(t2 + t4)
2392  end function sll_f_cubic_spline_2d_eval
2393 
2394  ! sll_f_cubic_spline_2d_eval_deriv_x1(): given discrete data f(i,j) that are
2395  ! described by a 2-dimensional cubic spline fit s(x1,x2), where the
2396  ! continuous variables x1 and x2 are within the original limits of i and j
2397  ! respectively, interpolate_x1_derivative() returns the value of
2398  !
2399  ! partial s
2400  ! -------------
2401  ! partial x1
2402  !
2403  ! evaluated at the point (x1,x2). (Sorry for the ambiguous use of x1)
2404 
2412  function sll_f_cubic_spline_2d_eval_deriv_x1(x1, x2, spline)
2414  intrinsic :: associated, int, real
2415  sll_real64, intent(in) :: x1
2416  sll_real64, intent(in) :: x2
2417  type(sll_t_cubic_spline_2d) :: spline
2418  sll_real64 :: rh1 ! reciprocal of cell spacing
2419  sll_real64 :: rh2 ! reciprocal of cell spacing
2420  sll_int32 :: cell
2421  sll_real64 :: dx
2422  sll_real64 :: cdx ! 1-dx
2423  sll_real64 :: t0 ! temp/scratch variables ...
2424  sll_real64 :: t1
2425  sll_real64 :: t2
2426  sll_real64 :: t3
2427  sll_real64 :: t4
2428  sll_real64 :: cim1 ! C_(i-1)
2429  sll_real64 :: ci ! C_i
2430  sll_real64 :: cip1 ! C_(i+1)
2431  sll_real64 :: cip2 ! C_(i+2)
2432  sll_real64 :: x1_min
2433  sll_real64 :: x2_min
2434  sll_int32 :: num_pts_x1
2435  sll_int32 :: num_pts_x2
2436  sll_real64, dimension(:), pointer :: coeffs_line_jm1
2437  sll_real64, dimension(:), pointer :: coeffs_line_j
2438  sll_real64, dimension(:), pointer :: coeffs_line_jp1
2439  sll_real64, dimension(:), pointer :: coeffs_line_jp2
2440  ! We set these as assertions since we want the flexibility of turning
2441  ! them off.
2442  sll_assert((x1 .ge. spline%x1_min) .and. (x1 .le. spline%x1_max))
2443  sll_assert((x2 .ge. spline%x2_min) .and. (x2 .le. spline%x2_max))
2444 
2445  x1_min = spline%x1_min
2446  x2_min = spline%x2_min
2447  num_pts_x1 = spline%num_pts_x1
2448  num_pts_x2 = spline%num_pts_x2
2449  rh1 = spline%x1_rdelta
2450  rh2 = spline%x2_rdelta
2451  ! find the cell and offset for x2
2452  t0 = (x2 - x2_min)*rh2
2453  cell = int(t0) + 1
2454  dx = t0 - real(cell - 1, f64)
2455  cdx = 1.0_f64 - dx
2456  ! write (*,'(a,i8, a, f20.12)') 'cell = ', cell, ', dx = ', dx
2457  ! interpolate the coefficients along the line of constant x1. These
2458  ! computations are independent from one another. A little problem is
2459  ! the redundancy in the computation of the cell and offset along each
2460  ! of the constant x2 lines, as this will be done by each call of
2461  ! interpolate_value_aux(). This suggests that the proper refactoring
2462  ! of this function would have the cell and offset as arguments.
2463  coeffs_line_jm1 => spline%coeffs(0:num_pts_x1 + 2, cell - 1)
2464  coeffs_line_j => spline%coeffs(0:num_pts_x1 + 2, cell)
2465  coeffs_line_jp1 => spline%coeffs(0:num_pts_x1 + 2, cell + 1)
2466  coeffs_line_jp2 => spline%coeffs(0:num_pts_x1 + 2, cell + 2)
2467  cim1 = interpolate_derivative_aux(x1, x1_min, rh1, coeffs_line_jm1)
2468  ci = interpolate_derivative_aux(x1, x1_min, rh1, coeffs_line_j)
2469  cip1 = interpolate_derivative_aux(x1, x1_min, rh1, coeffs_line_jp1)
2470  cip2 = interpolate_derivative_aux(x1, x1_min, rh1, coeffs_line_jp2)
2471  t1 = 3.0_f64*ci
2472  t3 = 3.0_f64*cip1
2473  t2 = cdx*(cdx*(cdx*(cim1 - t1) + t1) + t1) + ci
2474  t4 = dx*(dx*(dx*(cip2 - t3) + t3) + t3) + cip1
2475  sll_f_cubic_spline_2d_eval_deriv_x1 = (1.0_f64/6.0_f64)*(t2 + t4)
2477 
2478  ! sll_f_cubic_spline_2d_eval_deriv_x2(): given discrete data f(i,j) that are
2479  ! described by a 2-dimensional cubic spline fit s(x1,x2), where the
2480  ! continuous variables x1 and x2 are within the original limits of i and j
2481  ! respectively, interpolate_x1_derivative() returns the value of
2482  !
2483  ! partial s
2484  ! -------------
2485  ! partial x2
2486  !
2487  ! evaluated at the point (x1,x2). (Sorry for the ambiguous use of x1)
2488 
2498  function sll_f_cubic_spline_2d_eval_deriv_x2(x1, x2, spline)
2500  intrinsic :: associated, int, real
2501  sll_real64, intent(in) :: x1
2502  sll_real64, intent(in) :: x2
2503  type(sll_t_cubic_spline_2d) :: spline
2504  sll_real64 :: rh1 ! reciprocal of cell spacing
2505  sll_real64 :: rh2 ! reciprocal of cell spacing
2506  sll_int32 :: cell
2507  sll_real64 :: dx
2508  sll_real64 :: cdx ! 1-dx
2509  sll_real64 :: t0 ! temp/scratch variables ...
2510  sll_real64 :: t1
2511  sll_real64 :: t2
2512  sll_real64 :: t3
2513  sll_real64 :: cim1 ! C_(i-1)
2514  sll_real64 :: ci ! C_i
2515  sll_real64 :: cip1 ! C_(i+1)
2516  sll_real64 :: cip2 ! C_(i+2)
2517  sll_real64 :: x1_min
2518  sll_real64 :: x2_min
2519  sll_int32 :: num_pts_x1
2520  sll_int32 :: num_pts_x2
2521  sll_real64, dimension(:), pointer :: coeffs_line_jm1
2522  sll_real64, dimension(:), pointer :: coeffs_line_j
2523  sll_real64, dimension(:), pointer :: coeffs_line_jp1
2524  sll_real64, dimension(:), pointer :: coeffs_line_jp2
2525  ! We set these as assertions since we want the flexibility of turning
2526  ! them off.
2527  sll_assert((x1 .ge. spline%x1_min) .and. (x1 .le. spline%x1_max))
2528  sll_assert((x2 .ge. spline%x2_min) .and. (x2 .le. spline%x2_max))
2529 
2530  x1_min = spline%x1_min
2531  x2_min = spline%x2_min
2532  num_pts_x1 = spline%num_pts_x1
2533  num_pts_x2 = spline%num_pts_x2
2534  rh1 = spline%x1_rdelta
2535  rh2 = spline%x2_rdelta
2536  ! find the cell and offset for x2
2537  t0 = (x2 - x2_min)*rh2
2538  cell = int(t0) + 1
2539  dx = t0 - real(cell - 1, f64)
2540  cdx = 1.0_f64 - dx
2541  ! write (*,'(a,i8, a, f20.12)') 'cell = ', cell, ', dx = ', dx
2542  ! interpolate the coefficients along the line of constant x1. These
2543  ! computations are independent from one another. A little problem is
2544  ! the redundancy in the computation of the cell and offset along each
2545  ! of the constant x2 lines, as this will be done by each call of
2546  ! interpolate_value_aux(). This suggests that the proper refactoring
2547  ! of this function would have the cell and offset as arguments.
2548  coeffs_line_jm1 => spline%coeffs(0:num_pts_x1 + 2, cell - 1)
2549  coeffs_line_j => spline%coeffs(0:num_pts_x1 + 2, cell)
2550  coeffs_line_jp1 => spline%coeffs(0:num_pts_x1 + 2, cell + 1)
2551  coeffs_line_jp2 => spline%coeffs(0:num_pts_x1 + 2, cell + 2)
2552  cim1 = interpolate_value_aux(x1, x1_min, rh1, coeffs_line_jm1)
2553  ci = interpolate_value_aux(x1, x1_min, rh1, coeffs_line_j)
2554  cip1 = interpolate_value_aux(x1, x1_min, rh1, coeffs_line_jp1)
2555  cip2 = interpolate_value_aux(x1, x1_min, rh1, coeffs_line_jp2)
2556  t1 = 2.0_f64*(cim1 - 2.0_f64*ci + cip1)
2557  t2 = -cim1 + 3.0_f64*(ci - cip1) + cip2
2558  t3 = cip1 - cim1
2559  sll_f_cubic_spline_2d_eval_deriv_x2 = 0.5_f64*rh2*(dx*(t1 + dx*t2) + t3)
2561 
2562  subroutine sll_s_cubic_spline_2d_get_coeff(spline, coeff)
2563  type(sll_t_cubic_spline_2d) :: spline
2564  sll_real64, dimension(:), intent(out) :: coeff
2565  sll_int32 :: i
2566  sll_int32 :: j
2567  sll_int32 :: num_pts_x1
2568  sll_int32 :: num_pts_x2
2569 
2570  num_pts_x1 = spline%num_pts_x1
2571  num_pts_x2 = spline%num_pts_x2
2572 
2573  sll_assert(size(coeff) >= (num_pts_x1 + 2)*(num_pts_x2 + 2))
2574 
2575  do j = 1, num_pts_x2 + 2
2576  do i = 1, num_pts_x1 + 2
2577  coeff(i + (num_pts_x1 + 2)*(j - 1)) = spline%coeffs(i - 1, j - 1)
2578  end do
2579  end do
2580  end subroutine sll_s_cubic_spline_2d_get_coeff
2581 
2582  subroutine sll_s_cubic_spline_2d_free(spline)
2583  type(sll_t_cubic_spline_2d) :: spline
2584  sll_int32 :: ierr
2585  sll_deallocate(spline%d1, ierr)
2586  sll_deallocate(spline%d2, ierr)
2587  sll_deallocate(spline%coeffs, ierr)
2588  spline%data => null()
2589  if (associated(spline%x1_min_slopes)) then
2590  sll_deallocate(spline%x1_min_slopes, ierr)
2591  end if
2592  if (associated(spline%x1_max_slopes)) then
2593  sll_deallocate(spline%x1_max_slopes, ierr)
2594  end if
2595  if (associated(spline%x2_min_slopes)) then
2596  sll_deallocate(spline%x2_min_slopes, ierr)
2597  end if
2598  if (associated(spline%x2_max_slopes)) then
2599  sll_deallocate(spline%x2_max_slopes, ierr)
2600  end if
2601  if (associated(spline%x1_min_slopes_coeffs)) then
2602  sll_deallocate(spline%x1_min_slopes_coeffs, ierr)
2603  end if
2604  if (associated(spline%x1_max_slopes_coeffs)) then
2605  sll_deallocate(spline%x1_max_slopes_coeffs, ierr)
2606  end if
2607  if (associated(spline%x2_min_slopes_coeffs)) then
2608  sll_deallocate(spline%x2_min_slopes_coeffs, ierr)
2609  end if
2610  if (associated(spline%x2_max_slopes_coeffs)) then
2611  sll_deallocate(spline%x2_max_slopes_coeffs, ierr)
2612  end if
2613  end subroutine sll_s_cubic_spline_2d_free
2614 
2616  subroutine sll_s_cubic_spline_1d_eval_disp(spline, alpha, output_array)
2617  type(sll_t_cubic_spline_1d), intent(in) :: spline
2618  sll_real64, intent(in) :: alpha
2619  sll_real64, intent(out) :: output_array(:)
2620 
2621  !local variables
2622  sll_real64 :: alpha0
2623  sll_int32 :: dcell
2624  sll_real64 :: alpha1
2625  sll_int32 :: i
2626  sll_int32 :: cell
2627  sll_int32 :: num_cells
2628 
2629  ! Check that we are excepting the right number of output values
2630  sll_assert(size(output_array) == spline%n_points)
2631 
2632  alpha0 = alpha*spline%rdelta
2633  dcell = floor(alpha0)
2634  alpha1 = alpha0 - real(dcell, f64)
2635  select case (spline%bc_type)
2636  case (sll_p_periodic)
2637  num_cells = spline%n_points - 1
2638  case (sll_p_hermite)
2639  num_cells = spline%n_points
2640  end select
2641 
2642  do i = max(1, 1 - dcell), min(num_cells, num_cells - dcell)
2643  cell = i + dcell!modulo( i + dcell -1, this%spline%n_points-1) +1
2644  call spline_interpolate_from_interpolant_cell_dx(spline, cell, alpha1, output_array(i))
2645  end do
2646 
2647  if (spline%bc_type == sll_p_periodic) then
2648  do i = 1, max(1, 1 - dcell) - 1
2649  cell = modulo(i + dcell - 1, num_cells) + 1
2650  call spline_interpolate_from_interpolant_cell_dx(spline, cell, alpha1, output_array(i))
2651 
2652  end do
2653 
2654  do i = min(num_cells, num_cells - dcell) + 1, num_cells
2655  cell = modulo(i + dcell - 1, num_cells) + 1
2656  call spline_interpolate_from_interpolant_cell_dx(spline, cell, alpha1, output_array(i))
2657 
2658  end do
2659 
2660  ! First and last point equal for periodic
2661  output_array(spline%n_points) = output_array(1)
2662  else
2663  alpha1 = 0.0_f64
2664  if (dcell < 0) then
2665  cell = 1
2666  call spline_interpolate_from_interpolant_cell_dx(spline, cell, alpha1, output_array(1))
2667 
2668  do i = 2, -dcell
2669  output_array(i) = output_array(1)
2670  end do
2671  else
2672  cell = spline%n_points
2673  call spline_interpolate_from_interpolant_cell_dx(spline, cell, alpha1, output_array(spline%n_points))
2674 
2675  do i = num_cells - dcell + 1, spline%n_points - 1
2676  output_array(i) = output_array(spline%n_points)
2677  end do
2678 
2679  end if
2680  end if
2681 
2682  end subroutine sll_s_cubic_spline_1d_eval_disp
2683 
2685  subroutine spline_interpolate_from_interpolant_cell_dx(spline, cell, dx, out)
2686  type(sll_t_cubic_spline_1d), intent(in) :: spline
2687  sll_int32, intent(in) :: cell
2688  sll_real64, intent(in) :: dx
2689  sll_real64, intent(out) :: out
2690  sll_real64 :: cdx ! 1-dx
2691  sll_real64 :: cim1 ! C_(i-1)
2692  sll_real64 :: ci ! C_i
2693  sll_real64 :: cip1 ! C_(i+1)
2694  sll_real64 :: cip2 ! C_(i+2)
2695  sll_real64 :: t1
2696  sll_real64 :: t2
2697  sll_real64 :: t3
2698  sll_real64 :: t4
2699 
2700  cdx = 1.0_f64 - dx
2701  cim1 = spline%coeffs(cell - 1)
2702  ci = spline%coeffs(cell)
2703  cip1 = spline%coeffs(cell + 1)
2704  cip2 = spline%coeffs(cell + 2)
2705  t1 = 3.0_f64*ci
2706  t3 = 3.0_f64*cip1
2707  t2 = cdx*(cdx*(cdx*(cim1 - t1) + t1) + t1) + ci
2708  t4 = dx*(dx*(dx*(cip2 - t3) + t3) + t3) + cip1
2709  out = inv_6*(t2 + t4)
2710 
2712 
2713 #undef NUM_TERMS
2714 end module sll_m_cubic_splines
2715 
Provides capabilities for data and derivative interpolation with cubic B-splines and different bounda...
real(kind=f64), parameter inv_6
subroutine interpolate_pointer_values(ptr_in, ptr_out, n, spline)
returns the values of the images of a collection of abscissae,
subroutine, public sll_s_cubic_spline_2d_compute_interpolant(data, spline)
Computes the spline coefficients for the given data. The coeffcients are first computed in the second...
subroutine, public sll_s_cubic_spline_2d_deposit_value(x1, x2, spline, a_out)
Updated distribution function at time .
real(kind=f64) function, public sll_f_cubic_spline_1d_eval_deriv(x, spline)
returns the value of the derivative at the image of the abscissa 'x', using the spline coefficients s...
real(kind=f64) function, public sll_f_cubic_spline_2d_eval_deriv_x2(x1, x2, spline)
Returns the interpolated value of the derivative.
subroutine compute_spline_2d_hrmt_hrmt(data, spline)
subroutine, public sll_s_cubic_spline_1d_compute_interpolant(f, spline)
Computes the cubic spline coefficients corresponding to the 1D data array 'f'.
subroutine, public sll_s_cubic_spline_1d_init(self, num_points, xmin, xmax, bc_type, sl, sr, fast_algorithm)
Returns a pointer to a heap-allocated cubic spline object.
subroutine compute_spline_1d_hermite_aux(f, num_pts, d, slope_l, slope_r, delta, coeffs)
subroutine, public sll_s_cubic_spline_2d_get_coeff(spline, coeff)
subroutine, public sll_s_cubic_spline_1d_eval_disp(spline, alpha, output_array)
Computes the interpolated values at each grid point replaced by alpha for the precomputed spline coef...
subroutine interpolate_pointer_derivatives(ptr_in, ptr_out, num_pts, spline)
analogous to the sll_s_cubic_spline_1d_eval_array_deriv() subroutine but its input and output arrays ...
real(kind=f64) function interpolate_derivative_aux(x, xmin, rh, coeffs)
subroutine spline_interpolate_from_interpolant_cell_dx(spline, cell, dx, out)
Helper function for sll_s_cubic_spline_1d_eval_disp: evaluate spline in given cell and normalized dis...
subroutine, public sll_s_cubic_spline_1d_eval_array(a_in, a_out, n, spline)
returns the values of the images of a collection of abscissae, represented by a 1D array in another o...
subroutine compute_spline_2d_hrmt_prdc(data, spline)
subroutine compute_spline_1d_periodic(f, spline)
subroutine compute_spline_1d_hermite(f, spline)
subroutine, public sll_s_cubic_spline_1d_eval_array_deriv(array_in, array_out, num_pts, spline)
returns the values of the derivatives evaluated at a collection of abscissae stored in the input 1D a...
subroutine compute_spline_1d_periodic_aux(f, num_pts, d, coeffs)
subroutine compute_spline_2d_prdc_prdc(data, spline)
real(kind=f64) function, public sll_f_cubic_spline_1d_eval(x, spline)
returns the value of the interpolated image of the abscissa x, using the spline coefficients stored i...
real(kind=f64) function, public sll_f_cubic_spline_2d_eval_deriv_x1(x1, x2, spline)
Returns the interpolated value of the derivative in the x1 direction at the point (x1,...
real(kind=f64) function, public sll_f_cubic_spline_2d_eval(x1, x2, spline)
Returns the interpolated value of the image of the point (x1,x2) using the spline decomposition store...
real(kind=f64) function interpolate_value_aux(x, xmin, rh, coeffs)
subroutine, public sll_s_cubic_spline_1d_free(spline)
deallocate the sll_t_cubic_spline_1d object
subroutine, public sll_s_cubic_spline_2d_init(this, num_pts_x1, num_pts_x2, x1_min, x1_max, x2_min, x2_max, x1_bc_type, x2_bc_type, const_slope_x1_min, const_slope_x1_max, const_slope_x2_min, const_slope_x2_max, x1_min_slopes, x1_max_slopes, x2_min_slopes, x2_max_slopes)
subroutine compute_spline_2d_prdc_hrmt(data, spline)
type(sll_t_cubic_spline_2d) function, pointer sll_f_new_cubic_spline_2d(num_pts_x1, num_pts_x2, x1_min, x1_max, x2_min, x2_max, x1_bc_type, x2_bc_type, const_slope_x1_min, const_slope_x1_max, const_slope_x2_min, const_slope_x2_max, x1_min_slopes, x1_max_slopes, x2_min_slopes, x2_max_slopes)
Returns a pointer to a heap-allocated 2D cubic spline object.
subroutine, public sll_s_cubic_spline_2d_free(spline)
Tridiagonal system solver.
subroutine, public sll_s_setup_cyclic_tridiag(a, n, cts, ipiv)
Give the factorization of the matrix in argument.
subroutine, public sll_s_solve_cyclic_tridiag_double(cts, ipiv, b, n, x)
Solves tridiagonal system.
basic type for one-dimensional cubic spline data.
Basic type for one-dimensional cubic spline data.
    Report Typos and Errors