From 76563337055ed2814a0b36ade8b876b528faa458 Mon Sep 17 00:00:00 2001 From: Nick Wogan Date: Mon, 17 Jun 2024 12:34:10 -0700 Subject: [PATCH] adjusted linear interpolation module to get derivative --- .../linear_interpolation_module.F90 | 1302 ++--------------- 1 file changed, 120 insertions(+), 1182 deletions(-) diff --git a/src/dependencies/linear_interpolation_module.F90 b/src/dependencies/linear_interpolation_module.F90 index fabed69..554c55d 100755 --- a/src/dependencies/linear_interpolation_module.F90 +++ b/src/dependencies/linear_interpolation_module.F90 @@ -13,15 +13,10 @@ !@note The default real kind (`wp`) can be ! changed using optional preprocessor flags. ! This library was built with real kind: -#ifdef REAL32 -! `real(kind=real32)` [4 bytes] -#elif REAL64 -! `real(kind=real64)` [8 bytes] -#elif REAL128 -! `real(kind=real128)` [16 bytes] -#else -! `real(kind=real64)` [8 bytes] -#endif + +#:set TYPES = ['real(wp)', 'type(dual)'] +#:set NAMES = ['real', 'dual'] +#:set TYPES_NAMES = list(zip(TYPES, NAMES)) module linear_interpolation_module @@ -31,17 +26,7 @@ module linear_interpolation_module private -#ifdef REAL32 - integer,parameter,public :: finterp_rk = real32 !! real kind used by this module [4 bytes] -#elif REAL64 - integer,parameter,public :: finterp_rk = real64 !! real kind used by this module [8 bytes] -#elif REAL128 - integer,parameter,public :: finterp_rk = real128 !! real kind used by this module [16 bytes] -#else - integer,parameter,public :: finterp_rk = real64 !! real kind used by this module [8 bytes] -#endif - - integer,parameter :: wp = finterp_rk !! local copy of `finterp_rk` with a shorter name + integer,parameter :: wp = real64 !! local copy of `finterp_rk` with a shorter name real(wp),parameter,private :: zero = 0.0_wp !! numeric constant real(wp),parameter,private :: one = 1.0_wp !! numeric constant @@ -73,7 +58,8 @@ end subroutine destroy_func contains private procedure,public :: initialize => initialize_1d - procedure,public :: evaluate => interp_1d + procedure,public :: evaluate => interp_1d_real + procedure,public :: evaluate_derivative => interp_1d_derivative procedure,public :: destroy => destroy_1d final :: finalize_1d end type linear_interp_1d @@ -94,125 +80,9 @@ end subroutine destroy_func final :: finalize_2d end type linear_interp_2d - type,extends(linear_interp_class),public :: linear_interp_3d - !! Class for 3d linear interpolation. - private - real(wp),dimension(:,:,:),allocatable :: f - real(wp),dimension(:),allocatable :: x - real(wp),dimension(:),allocatable :: y - real(wp),dimension(:),allocatable :: z - integer :: ilox = 1 - integer :: iloy = 1 - integer :: iloz = 1 - contains - private - procedure,public :: initialize => initialize_3d - procedure,public :: evaluate => interp_3d - procedure,public :: destroy => destroy_3d - final :: finalize_3d - end type linear_interp_3d - - type,extends(linear_interp_class),public :: linear_interp_4d - !! Class for 4d linear interpolation. - private - real(wp),dimension(:,:,:,:),allocatable :: f - real(wp),dimension(:),allocatable :: x - real(wp),dimension(:),allocatable :: y - real(wp),dimension(:),allocatable :: z - real(wp),dimension(:),allocatable :: q - integer :: ilox = 1 - integer :: iloy = 1 - integer :: iloz = 1 - integer :: iloq = 1 - contains - private - procedure,public :: initialize => initialize_4d - procedure,public :: evaluate => interp_4d - procedure,public :: destroy => destroy_4d - final :: finalize_4d - end type linear_interp_4d - - type,extends(linear_interp_class),public :: linear_interp_5d - !! Class for 5d linear interpolation. - private - real(wp),dimension(:,:,:,:,:),allocatable :: f - real(wp),dimension(:),allocatable :: x - real(wp),dimension(:),allocatable :: y - real(wp),dimension(:),allocatable :: z - real(wp),dimension(:),allocatable :: q - real(wp),dimension(:),allocatable :: r - integer :: ilox = 1 - integer :: iloy = 1 - integer :: iloz = 1 - integer :: iloq = 1 - integer :: ilor = 1 - contains - private - procedure,public :: initialize => initialize_5d - procedure,public :: evaluate => interp_5d - procedure,public :: destroy => destroy_5d - final :: finalize_5d - end type linear_interp_5d - - type,extends(linear_interp_class),public :: linear_interp_6d - !! Class for 6d linear interpolation. - private - real(wp),dimension(:,:,:,:,:,:),allocatable :: f - real(wp),dimension(:),allocatable :: x - real(wp),dimension(:),allocatable :: y - real(wp),dimension(:),allocatable :: z - real(wp),dimension(:),allocatable :: q - real(wp),dimension(:),allocatable :: r - real(wp),dimension(:),allocatable :: s - integer :: ilox = 1 - integer :: iloy = 1 - integer :: iloz = 1 - integer :: iloq = 1 - integer :: ilor = 1 - integer :: ilos = 1 - contains - private - procedure,public :: initialize => initialize_6d - procedure,public :: evaluate => interp_6d - procedure,public :: destroy => destroy_6d - final :: finalize_6d - end type linear_interp_6d - - type,extends(linear_interp_1d),public :: nearest_interp_1d - !! Class for 1d nearest neighbor interpolation. - contains - procedure,public :: evaluate => nearest_1d - end type nearest_interp_1d - - type,extends(linear_interp_2d),public :: nearest_interp_2d - !! Class for 2d nearest neighbor interpolation. - contains - procedure,public :: evaluate => nearest_2d - end type nearest_interp_2d - - type,extends(linear_interp_3d),public :: nearest_interp_3d - !! Class for 3d nearest neighbor interpolation. - contains - procedure,public :: evaluate => nearest_3d - end type nearest_interp_3d - - type,extends(linear_interp_4d),public :: nearest_interp_4d - !! Class for 4d nearest neighbor interpolation. - contains - procedure,public :: evaluate => nearest_4d - end type nearest_interp_4d - - type,extends(linear_interp_5d),public :: nearest_interp_5d - !! Class for 5d nearest neighbor interpolation. - contains - procedure,public :: evaluate => nearest_5d - end type nearest_interp_5d - - type,extends(linear_interp_6d),public :: nearest_interp_6d - !! Class for 6d nearest neighbor interpolation. - contains - procedure,public :: evaluate => nearest_6d - end type nearest_interp_6d + interface dintrv + module procedure :: dintrv_real, dintrv_dual + end interface contains !***************************************************************************************** @@ -245,62 +115,6 @@ pure elemental subroutine finalize_2d(me) end subroutine finalize_2d !***************************************************************************************** -!***************************************************************************************** -!> -! Finalizer for a [[linear_interp_3d]] type. - - pure elemental subroutine finalize_3d(me) - - implicit none - - type(linear_interp_3d),intent(inout) :: me - call me%destroy() - - end subroutine finalize_3d -!***************************************************************************************** - -!***************************************************************************************** -!> -! Finalizer for a [[linear_interp_4d]] type. - - pure elemental subroutine finalize_4d(me) - - implicit none - - type(linear_interp_4d),intent(inout) :: me - call me%destroy() - - end subroutine finalize_4d -!***************************************************************************************** - -!***************************************************************************************** -!> -! Finalizer for a [[linear_interp_5d]] type. - - pure elemental subroutine finalize_5d(me) - - implicit none - - type(linear_interp_5d),intent(inout) :: me - call me%destroy() - - end subroutine finalize_5d -!***************************************************************************************** - -!***************************************************************************************** -!> -! Finalizer for a [[linear_interp_6d]] type. - - pure elemental subroutine finalize_6d(me) - - implicit none - - type(linear_interp_6d),intent(inout) :: me - call me%destroy() - - end subroutine finalize_6d -!***************************************************************************************** - !***************************************************************************************** !> ! Destructor for a [[linear_interp_1d]] class. @@ -341,810 +155,176 @@ end subroutine destroy_2d !***************************************************************************************** !> -! Destructor for a [[linear_interp_3d]] class. - - pure elemental subroutine destroy_3d(me) - - implicit none - - class(linear_interp_3d),intent(inout) :: me - - if (allocated(me%f)) deallocate(me%f) - if (allocated(me%x)) deallocate(me%x) - if (allocated(me%y)) deallocate(me%y) - if (allocated(me%z)) deallocate(me%z) - me%ilox = 1 - me%iloy = 1 - me%iloz = 1 - me%initialized = .false. - - end subroutine destroy_3d -!***************************************************************************************** - -!***************************************************************************************** -!> -! Destructor for a [[linear_interp_4d]] class. - - pure elemental subroutine destroy_4d(me) - - implicit none - - class(linear_interp_4d),intent(inout) :: me - - if (allocated(me%f)) deallocate(me%f) - if (allocated(me%x)) deallocate(me%x) - if (allocated(me%y)) deallocate(me%y) - if (allocated(me%z)) deallocate(me%z) - if (allocated(me%q)) deallocate(me%q) - me%ilox = 1 - me%iloy = 1 - me%iloz = 1 - me%iloq = 1 - me%initialized = .false. - - end subroutine destroy_4d -!***************************************************************************************** - -!***************************************************************************************** -!> -! Destructor for a [[linear_interp_5d]] class. - - pure elemental subroutine destroy_5d(me) - - implicit none - - class(linear_interp_5d),intent(inout) :: me - - if (allocated(me%f)) deallocate(me%f) - if (allocated(me%x)) deallocate(me%x) - if (allocated(me%y)) deallocate(me%y) - if (allocated(me%z)) deallocate(me%z) - if (allocated(me%q)) deallocate(me%q) - if (allocated(me%r)) deallocate(me%r) - me%ilox = 1 - me%iloy = 1 - me%iloz = 1 - me%iloq = 1 - me%ilor = 1 - me%initialized = .false. - - end subroutine destroy_5d -!***************************************************************************************** - -!***************************************************************************************** -!> -! Destructor for a [[linear_interp_6d]] class. - - pure elemental subroutine destroy_6d(me) - - implicit none - - class(linear_interp_6d),intent(inout) :: me - - if (allocated(me%f)) deallocate(me%f) - if (allocated(me%x)) deallocate(me%x) - if (allocated(me%y)) deallocate(me%y) - if (allocated(me%z)) deallocate(me%z) - if (allocated(me%q)) deallocate(me%q) - if (allocated(me%r)) deallocate(me%r) - if (allocated(me%s)) deallocate(me%s) - me%ilox = 1 - me%iloy = 1 - me%iloz = 1 - me%iloq = 1 - me%ilor = 1 - me%ilos = 1 - me%initialized = .false. - - end subroutine destroy_6d -!***************************************************************************************** - -!***************************************************************************************** -!> -! Constructor for a [[linear_interp_1d]] class. - - pure subroutine initialize_1d(me,x,f,istat) - - implicit none - - class(linear_interp_1d),intent(inout) :: me - real(wp),dimension(:),intent(in) :: x - real(wp),dimension(:),intent(in) :: f - integer,intent(out) :: istat !! `0` : no problems, - !! `1` : `x` is not strictly increasing, - !! `10` : `x` is not equal to size(f,1), - !! `100` : cannot use linear interpolation for only one point. - - call me%destroy() - - istat = 0 - - if (istat==0 .and. size(x)/=size(f,1)) istat = 10 - - if (istat==0) then - call me%check_inputs(x=x,ierr=istat) - if (istat==0) then - allocate(me%f(size(x))); me%f = f - allocate(me%x(size(x))); me%x = x - me%initialized = .true. - end if - end if - - end subroutine initialize_1d -!***************************************************************************************** - -!***************************************************************************************** -!> -! Constructor for a [[linear_interp_2d]] class. - - pure subroutine initialize_2d(me,x,y,f,istat) - - implicit none - - class(linear_interp_2d),intent(inout) :: me - real(wp),dimension(:),intent(in) :: x - real(wp),dimension(:),intent(in) :: y - real(wp),dimension(:,:),intent(in) :: f - integer,intent(out) :: istat !! `0` : no problems, - !! `1` : `x` is not strictly increasing, - !! `2` : `y` is not strictly increasing, - !! `10` : `x` is not equal to size(f,1), - !! `20` : `y` is not equal to size(f,2), - !! `100` : cannot use linear interpolation for only one point. - - call me%destroy() - - istat = 0 - - if (istat==0 .and. size(x)/=size(f,1)) istat = 10 - if (istat==0 .and. size(y)/=size(f,2)) istat = 20 - - if (istat==0) then - call me%check_inputs(x=x,y=y,ierr=istat) - if (istat==0) then - allocate(me%f(size(x),size(y))); me%f = f - allocate(me%x(size(x))); me%x = x - allocate(me%y(size(y))); me%y = y - me%initialized = .true. - end if - end if - - end subroutine initialize_2d -!***************************************************************************************** - -!***************************************************************************************** -!> -! Constructor for a [[linear_interp_3d]] class. - - pure subroutine initialize_3d(me,x,y,z,f,istat) - - implicit none - - class(linear_interp_3d),intent(inout) :: me - real(wp),dimension(:),intent(in) :: x - real(wp),dimension(:),intent(in) :: y - real(wp),dimension(:),intent(in) :: z - real(wp),dimension(:,:,:),intent(in) :: f - integer,intent(out) :: istat !! `0` : no problems, - !! `1` : `x` is not strictly increasing, - !! `2` : `y` is not strictly increasing, - !! `3` : `z` is not strictly increasing, - !! `10` : `x` is not equal to size(f,1), - !! `20` : `y` is not equal to size(f,2), - !! `30` : `z` is not equal to size(f,3), - !! `100` : cannot use linear interpolation for only one point. - - call me%destroy() - - istat = 0 - - if (istat==0 .and. size(x)/=size(f,1)) istat = 10 - if (istat==0 .and. size(y)/=size(f,2)) istat = 20 - if (istat==0 .and. size(z)/=size(f,3)) istat = 30 - - if (istat==0) then - call me%check_inputs(x=x,y=y,z=z,ierr=istat) - if (istat==0) then - allocate(me%f(size(x),size(y),size(z))); me%f = f - allocate(me%x(size(x))); me%x = x - allocate(me%y(size(y))); me%y = y - allocate(me%z(size(z))); me%z = z - me%initialized = .true. - end if - end if - - end subroutine initialize_3d -!***************************************************************************************** - -!***************************************************************************************** -!> -! Constructor for a [[linear_interp_4d]] class. - - pure subroutine initialize_4d(me,x,y,z,q,f,istat) - - implicit none - - class(linear_interp_4d),intent(inout) :: me - real(wp),dimension(:),intent(in) :: x - real(wp),dimension(:),intent(in) :: y - real(wp),dimension(:),intent(in) :: z - real(wp),dimension(:),intent(in) :: q - real(wp),dimension(:,:,:,:),intent(in) :: f - integer,intent(out) :: istat !! `0` : no problems, - !! `1` : `x` is not strictly increasing, - !! `2` : `y` is not strictly increasing, - !! `3` : `z` is not strictly increasing, - !! `4` : `q` is not strictly increasing, - !! `10` : `x` is not equal to size(f,1), - !! `20` : `y` is not equal to size(f,2), - !! `30` : `z` is not equal to size(f,3), - !! `40` : `q` is not equal to size(f,4), - !! `100` : cannot use linear interpolation for only one point. - - call me%destroy() - - istat = 0 - - if (istat==0 .and. size(x)/=size(f,1)) istat = 10 - if (istat==0 .and. size(y)/=size(f,2)) istat = 20 - if (istat==0 .and. size(z)/=size(f,3)) istat = 30 - if (istat==0 .and. size(q)/=size(f,4)) istat = 40 - - if (istat==0) then - call me%check_inputs(x=x,y=y,z=z,q=q,ierr=istat) - if (istat==0) then - allocate(me%f(size(x),size(y),size(z),size(q))); me%f = f - allocate(me%x(size(x))); me%x = x - allocate(me%y(size(y))); me%y = y - allocate(me%z(size(z))); me%z = z - allocate(me%q(size(q))); me%q = q - me%initialized = .true. - end if - end if - - end subroutine initialize_4d -!***************************************************************************************** - -!***************************************************************************************** -!> -! Constructor for a [[linear_interp_5d]] class. - - pure subroutine initialize_5d(me,x,y,z,q,r,f,istat) - - implicit none - - class(linear_interp_5d),intent(inout) :: me - real(wp),dimension(:),intent(in) :: x - real(wp),dimension(:),intent(in) :: y - real(wp),dimension(:),intent(in) :: z - real(wp),dimension(:),intent(in) :: q - real(wp),dimension(:),intent(in) :: r - real(wp),dimension(:,:,:,:,:),intent(in) :: f - integer,intent(out) :: istat !! `0` : no problems, - !! `1` : `x` is not strictly increasing, - !! `2` : `y` is not strictly increasing, - !! `3` : `z` is not strictly increasing, - !! `4` : `q` is not strictly increasing, - !! `5` : `r` is not strictly increasing, - !! `10` : `x` is not equal to size(f,1), - !! `20` : `y` is not equal to size(f,2), - !! `30` : `z` is not equal to size(f,3), - !! `40` : `q` is not equal to size(f,4), - !! `50` : `r` is not equal to size(f,5), - !! `100` : cannot use linear interpolation for only one point. - - call me%destroy() - - istat = 0 - - if (istat==0 .and. size(x)/=size(f,1)) istat = 10 - if (istat==0 .and. size(y)/=size(f,2)) istat = 20 - if (istat==0 .and. size(z)/=size(f,3)) istat = 30 - if (istat==0 .and. size(q)/=size(f,4)) istat = 40 - if (istat==0 .and. size(r)/=size(f,5)) istat = 50 - - if (istat==0) then - call me%check_inputs(x=x,y=y,z=z,q=q,r=r,ierr=istat) - if (istat==0) then - allocate(me%f(size(x),size(y),size(z),size(q),size(r))); me%f = f - allocate(me%x(size(x))); me%x = x - allocate(me%y(size(y))); me%y = y - allocate(me%z(size(z))); me%z = z - allocate(me%q(size(q))); me%q = q - allocate(me%r(size(r))); me%r = r - me%initialized = .true. - end if - end if - - end subroutine initialize_5d -!***************************************************************************************** - -!***************************************************************************************** -!> -! Constructor for a [[linear_interp_6d]] class. - - pure subroutine initialize_6d(me,x,y,z,q,r,s,f,istat) - - implicit none - - class(linear_interp_6d),intent(inout) :: me - real(wp),dimension(:),intent(in) :: x - real(wp),dimension(:),intent(in) :: y - real(wp),dimension(:),intent(in) :: z - real(wp),dimension(:),intent(in) :: q - real(wp),dimension(:),intent(in) :: r - real(wp),dimension(:),intent(in) :: s - real(wp),dimension(:,:,:,:,:,:),intent(in) :: f - integer,intent(out) :: istat !! `0` : no problems, - !! `1` : `x` is not strictly increasing, - !! `2` : `y` is not strictly increasing, - !! `3` : `z` is not strictly increasing, - !! `4` : `q` is not strictly increasing, - !! `5` : `r` is not strictly increasing, - !! `6` : `s` is not strictly increasing, - !! `10` : `x` is not equal to size(f,1), - !! `20` : `y` is not equal to size(f,2), - !! `30` : `z` is not equal to size(f,3), - !! `40` : `q` is not equal to size(f,4), - !! `50` : `r` is not equal to size(f,5), - !! `60` : `s` is not equal to size(f,6), - !! `100` : cannot use linear interpolation for only one point. - - call me%destroy() - - istat = 0 - - if (istat==0 .and. size(x)/=size(f,1)) istat = 10 - if (istat==0 .and. size(y)/=size(f,2)) istat = 20 - if (istat==0 .and. size(z)/=size(f,3)) istat = 30 - if (istat==0 .and. size(q)/=size(f,4)) istat = 40 - if (istat==0 .and. size(r)/=size(f,5)) istat = 50 - if (istat==0 .and. size(s)/=size(f,6)) istat = 60 - - if (istat==0) then - call me%check_inputs(x=x,y=y,z=z,q=q,r=r,s=s,ierr=istat) - if (istat==0) then - allocate(me%f(size(x),size(y),size(z),size(q),size(r),size(s))); me%f = f - allocate(me%x(size(x))); me%x = x - allocate(me%y(size(y))); me%y = y - allocate(me%z(size(z))); me%z = z - allocate(me%q(size(q))); me%q = q - allocate(me%r(size(r))); me%r = r - allocate(me%s(size(s))); me%s = s - me%initialized = .true. - end if - end if - - end subroutine initialize_6d -!***************************************************************************************** - -!***************************************************************************************** -!> -! 1D linear interpolation routine. - - pure subroutine interp_1d(me,x,f,istat) - - implicit none - - class(linear_interp_1d),intent(inout) :: me - real(wp),intent(in) :: x - real(wp),intent(out) :: f !! Interpolated \( f(x) \) - integer,intent(out),optional :: istat !! `0` : no problems, - !! `-1` : class has not been initialized - - integer,dimension(2) :: ix - real(wp) :: p1 - real(wp) :: q1 - integer :: mflag - - if (me%initialized) then - - call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag) - - q1 = (x-me%x(ix(1)))/(me%x(ix(2))-me%x(ix(1))) - p1 = one-q1 - - f = p1*me%f(ix(1)) + q1*me%f(ix(2)) - if (present(istat)) istat = 0 - - else - - if (present(istat)) istat = -1 - f = zero - - end if - - end subroutine interp_1d -!***************************************************************************************** - -!***************************************************************************************** -!> -! 2D linear interpolation routine. - - pure subroutine interp_2d(me,x,y,f,istat) - - implicit none - - class(linear_interp_2d),intent(inout) :: me - real(wp),intent(in) :: x - real(wp),intent(in) :: y - real(wp),intent(out) :: f !! Interpolated \( f(x,y) \) - integer,intent(out),optional :: istat !! `0` : no problems, - !! `-1` : class has not been initialized - - integer,dimension(2) :: ix, iy - real(wp) :: p1, p2 - real(wp) :: q1, q2 - integer :: mflag - real(wp) :: fx1, fx2 - - if (me%initialized) then - - call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag) - call dintrv(me%y,y,me%iloy,iy(1),iy(2),mflag) - - q1 = (x-me%x(ix(1)))/(me%x(ix(2))-me%x(ix(1))) - q2 = (y-me%y(iy(1)))/(me%y(iy(2))-me%y(iy(1))) - p1 = one-q1 - p2 = one-q2 - - fx1 = p1*me%f(ix(1),iy(1)) + q1*me%f(ix(2),iy(1)) - fx2 = p1*me%f(ix(1),iy(2)) + q1*me%f(ix(2),iy(2)) - - f = p2*fx1 + q2*fx2 - if (present(istat)) istat = 0 - - else - - if (present(istat)) istat = -1 - f = zero - - end if - - end subroutine interp_2d -!***************************************************************************************** - -!***************************************************************************************** -!> -! 3D linear interpolation routine. +! Constructor for a [[linear_interp_1d]] class. - pure subroutine interp_3d(me,x,y,z,f,istat) + pure subroutine initialize_1d(me,x,f,istat) implicit none - class(linear_interp_3d),intent(inout) :: me - real(wp),intent(in) :: x - real(wp),intent(in) :: y - real(wp),intent(in) :: z - real(wp),intent(out) :: f !! Interpolated \( f(x,y,z) \) - integer,intent(out),optional :: istat !! `0` : no problems, - !! `-1` : class has not been initialized - - integer,dimension(2) :: ix, iy, iz - real(wp) :: p1, p2, p3 - real(wp) :: q1, q2, q3 - integer :: mflag - real(wp) :: fx11, fx21, fx12, fx22, fxy1, fxy2 - - if (me%initialized) then - - call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag) - call dintrv(me%y,y,me%iloy,iy(1),iy(2),mflag) - call dintrv(me%z,z,me%iloz,iz(1),iz(2),mflag) - - q1 = (x-me%x(ix(1)))/(me%x(ix(2))-me%x(ix(1))) - q2 = (y-me%y(iy(1)))/(me%y(iy(2))-me%y(iy(1))) - q3 = (z-me%z(iz(1)))/(me%z(iz(2))-me%z(iz(1))) - p1 = one-q1 - p2 = one-q2 - p3 = one-q3 - - fx11 = p1*me%f(ix(1),iy(1),iz(1)) + q1*me%f(ix(2),iy(1),iz(1)) - fx21 = p1*me%f(ix(1),iy(2),iz(1)) + q1*me%f(ix(2),iy(2),iz(1)) - fx12 = p1*me%f(ix(1),iy(1),iz(2)) + q1*me%f(ix(2),iy(1),iz(2)) - fx22 = p1*me%f(ix(1),iy(2),iz(2)) + q1*me%f(ix(2),iy(2),iz(2)) - fxy1 = p2*fx11 + q2*fx21 - fxy2 = p2*fx12 + q2*fx22 + class(linear_interp_1d),intent(inout) :: me + real(wp),dimension(:),intent(in) :: x + real(wp),dimension(:),intent(in) :: f + integer,intent(out) :: istat !! `0` : no problems, + !! `1` : `x` is not strictly increasing, + !! `10` : `x` is not equal to size(f,1), + !! `100` : cannot use linear interpolation for only one point. - f = p3*fxy1 + q3*fxy2 - if (present(istat)) istat = 0 + call me%destroy() - else + istat = 0 - if (present(istat)) istat = -1 - f = zero + if (istat==0 .and. size(x)/=size(f,1)) istat = 10 + if (istat==0) then + call me%check_inputs(x=x,ierr=istat) + if (istat==0) then + allocate(me%f(size(x))); me%f = f + allocate(me%x(size(x))); me%x = x + me%initialized = .true. + end if end if - end subroutine interp_3d + end subroutine initialize_1d !***************************************************************************************** !***************************************************************************************** !> -! 4D linear interpolation routine. +! Constructor for a [[linear_interp_2d]] class. - pure subroutine interp_4d(me,x,y,z,q,f,istat) + pure subroutine initialize_2d(me,x,y,f,istat) implicit none - class(linear_interp_4d),intent(inout) :: me - real(wp),intent(in) :: x - real(wp),intent(in) :: y - real(wp),intent(in) :: z - real(wp),intent(in) :: q - real(wp),intent(out) :: f !! Interpolated \( f(x,y,z,q) \) - integer,intent(out),optional :: istat !! `0` : no problems, - !! `-1` : class has not been initialized - - integer,dimension(2) :: ix, iy, iz, iq - real(wp) :: p1, p2, p3, p4 - real(wp) :: q1, q2, q3, q4 - integer :: mflag - real(wp) :: fx111,fx211,fx121,fx221,fxy11,fxy21,fxyz1,& - fx112,fx212,fx122,fx222,fxy12,fxy22,fxyz2 - - if (me%initialized) then - - call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag) - call dintrv(me%y,y,me%iloy,iy(1),iy(2),mflag) - call dintrv(me%z,z,me%iloz,iz(1),iz(2),mflag) - call dintrv(me%q,q,me%iloq,iq(1),iq(2),mflag) + class(linear_interp_2d),intent(inout) :: me + real(wp),dimension(:),intent(in) :: x + real(wp),dimension(:),intent(in) :: y + real(wp),dimension(:,:),intent(in) :: f + integer,intent(out) :: istat !! `0` : no problems, + !! `1` : `x` is not strictly increasing, + !! `2` : `y` is not strictly increasing, + !! `10` : `x` is not equal to size(f,1), + !! `20` : `y` is not equal to size(f,2), + !! `100` : cannot use linear interpolation for only one point. - q1 = (x-me%x(ix(1)))/(me%x(ix(2))-me%x(ix(1))) - q2 = (y-me%y(iy(1)))/(me%y(iy(2))-me%y(iy(1))) - q3 = (z-me%z(iz(1)))/(me%z(iz(2))-me%z(iz(1))) - q4 = (q-me%q(iq(1)))/(me%q(iq(2))-me%q(iq(1))) - p1 = one-q1 - p2 = one-q2 - p3 = one-q3 - p4 = one-q4 - - fx111 = p1*me%f(ix(1),iy(1),iz(1),iq(1)) + q1*me%f(ix(2),iy(1),iz(1),iq(1)) - fx211 = p1*me%f(ix(1),iy(2),iz(1),iq(1)) + q1*me%f(ix(2),iy(2),iz(1),iq(1)) - fx121 = p1*me%f(ix(1),iy(1),iz(2),iq(1)) + q1*me%f(ix(2),iy(1),iz(2),iq(1)) - fx221 = p1*me%f(ix(1),iy(2),iz(2),iq(1)) + q1*me%f(ix(2),iy(2),iz(2),iq(1)) - fx112 = p1*me%f(ix(1),iy(1),iz(1),iq(2)) + q1*me%f(ix(2),iy(1),iz(1),iq(2)) - fx212 = p1*me%f(ix(1),iy(2),iz(1),iq(2)) + q1*me%f(ix(2),iy(2),iz(1),iq(2)) - fx122 = p1*me%f(ix(1),iy(1),iz(2),iq(2)) + q1*me%f(ix(2),iy(1),iz(2),iq(2)) - fx222 = p1*me%f(ix(1),iy(2),iz(2),iq(2)) + q1*me%f(ix(2),iy(2),iz(2),iq(2)) - - fxy11 = p2*fx111 + q2*fx211 - fxy21 = p2*fx121 + q2*fx221 - fxy12 = p2*fx112 + q2*fx212 - fxy22 = p2*fx122 + q2*fx222 - - fxyz1 = p3*fxy11 + q3*fxy21 - fxyz2 = p3*fxy12 + q3*fxy22 - - f = p4*fxyz1 + q4*fxyz2 - if (present(istat)) istat = 0 + call me%destroy() - else + istat = 0 - if (present(istat)) istat = -1 - f = zero + if (istat==0 .and. size(x)/=size(f,1)) istat = 10 + if (istat==0 .and. size(y)/=size(f,2)) istat = 20 + if (istat==0) then + call me%check_inputs(x=x,y=y,ierr=istat) + if (istat==0) then + allocate(me%f(size(x),size(y))); me%f = f + allocate(me%x(size(x))); me%x = x + allocate(me%y(size(y))); me%y = y + me%initialized = .true. + end if end if - end subroutine interp_4d + end subroutine initialize_2d !***************************************************************************************** !***************************************************************************************** !> -! 5D linear interpolation routine. - - pure subroutine interp_5d(me,x,y,z,q,r,f,istat) +! 1D linear interpolation routine. + #:for TYPE1, NAME in TYPES_NAMES + pure subroutine interp_1d_${NAME}$(me,x,f,istat) + #:if NAME == 'dual' + use forwarddiff + #:endif implicit none - class(linear_interp_5d),intent(inout) :: me - real(wp),intent(in) :: x - real(wp),intent(in) :: y - real(wp),intent(in) :: z - real(wp),intent(in) :: q - real(wp),intent(in) :: r - real(wp),intent(out) :: f !! Interpolated \( f(x,y,z,q,r) \) - integer,intent(out),optional :: istat !! `0` : no problems, - !! `-1` : class has not been initialized - - integer,dimension(2) :: ix, iy, iz, iq, ir - real(wp) :: p1, p2, p3, p4, p5 - real(wp) :: q1, q2, q3, q4, q5 + class(linear_interp_1d),intent(inout) :: me + ${TYPE1}$,intent(in) :: x + ${TYPE1}$,intent(out) :: f !! Interpolated \( f(x) \) + integer,intent(out),optional :: istat !! `0` : no problems, + !! `-1` : class has not been initialized + + integer,dimension(2) :: ix + ${TYPE1}$ :: p1 + ${TYPE1}$ :: q1 integer :: mflag - real(wp) :: fx1111, fx2111, fx1211, fx2211, fx1121, fx2121, fx1221, fx2221, & - fxy111, fxy211, fxy121, fxy221, fxyz11, fxyz21, fxyzq1, fx1112, & - fx2112, fx1212, fx2212, fx1122, fx2122, fx1222, fx2222, fxy112, & - fxy212, fxy122, fxy222, fxyz12, fxyz22, fxyzq2 if (me%initialized) then call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag) - call dintrv(me%y,y,me%iloy,iy(1),iy(2),mflag) - call dintrv(me%z,z,me%iloz,iz(1),iz(2),mflag) - call dintrv(me%q,q,me%iloq,iq(1),iq(2),mflag) - call dintrv(me%r,r,me%ilor,ir(1),ir(2),mflag) q1 = (x-me%x(ix(1)))/(me%x(ix(2))-me%x(ix(1))) - q2 = (y-me%y(iy(1)))/(me%y(iy(2))-me%y(iy(1))) - q3 = (z-me%z(iz(1)))/(me%z(iz(2))-me%z(iz(1))) - q4 = (q-me%q(iq(1)))/(me%q(iq(2))-me%q(iq(1))) - q5 = (r-me%r(ir(1)))/(me%r(ir(2))-me%r(ir(1))) p1 = one-q1 - p2 = one-q2 - p3 = one-q3 - p4 = one-q4 - p5 = one-q5 - - fx1111 = p1*me%f(ix(1),iy(1),iz(1),iq(1),ir(1)) + q1*me%f(ix(2),iy(1),iz(1),iq(1),ir(1)) - fx2111 = p1*me%f(ix(1),iy(2),iz(1),iq(1),ir(1)) + q1*me%f(ix(2),iy(2),iz(1),iq(1),ir(1)) - fx1211 = p1*me%f(ix(1),iy(1),iz(2),iq(1),ir(1)) + q1*me%f(ix(2),iy(1),iz(2),iq(1),ir(1)) - fx2211 = p1*me%f(ix(1),iy(2),iz(2),iq(1),ir(1)) + q1*me%f(ix(2),iy(2),iz(2),iq(1),ir(1)) - fx1121 = p1*me%f(ix(1),iy(1),iz(1),iq(2),ir(1)) + q1*me%f(ix(2),iy(1),iz(1),iq(2),ir(1)) - fx2121 = p1*me%f(ix(1),iy(2),iz(1),iq(2),ir(1)) + q1*me%f(ix(2),iy(2),iz(1),iq(2),ir(1)) - fx1221 = p1*me%f(ix(1),iy(1),iz(2),iq(2),ir(1)) + q1*me%f(ix(2),iy(1),iz(2),iq(2),ir(1)) - fx2221 = p1*me%f(ix(1),iy(2),iz(2),iq(2),ir(1)) + q1*me%f(ix(2),iy(2),iz(2),iq(2),ir(1)) - fx1112 = p1*me%f(ix(1),iy(1),iz(1),iq(1),ir(2)) + q1*me%f(ix(2),iy(1),iz(1),iq(1),ir(2)) - fx2112 = p1*me%f(ix(1),iy(2),iz(1),iq(1),ir(2)) + q1*me%f(ix(2),iy(2),iz(1),iq(1),ir(2)) - fx1212 = p1*me%f(ix(1),iy(1),iz(2),iq(1),ir(2)) + q1*me%f(ix(2),iy(1),iz(2),iq(1),ir(2)) - fx2212 = p1*me%f(ix(1),iy(2),iz(2),iq(1),ir(2)) + q1*me%f(ix(2),iy(2),iz(2),iq(1),ir(2)) - fx1122 = p1*me%f(ix(1),iy(1),iz(1),iq(2),ir(2)) + q1*me%f(ix(2),iy(1),iz(1),iq(2),ir(2)) - fx2122 = p1*me%f(ix(1),iy(2),iz(1),iq(2),ir(2)) + q1*me%f(ix(2),iy(2),iz(1),iq(2),ir(2)) - fx1222 = p1*me%f(ix(1),iy(1),iz(2),iq(2),ir(2)) + q1*me%f(ix(2),iy(1),iz(2),iq(2),ir(2)) - fx2222 = p1*me%f(ix(1),iy(2),iz(2),iq(2),ir(2)) + q1*me%f(ix(2),iy(2),iz(2),iq(2),ir(2)) - - fxy111 = p2*fx1111 + q2*fx2111 - fxy211 = p2*fx1211 + q2*fx2211 - fxy121 = p2*fx1121 + q2*fx2121 - fxy221 = p2*fx1221 + q2*fx2221 - fxy112 = p2*fx1112 + q2*fx2112 - fxy212 = p2*fx1212 + q2*fx2212 - fxy122 = p2*fx1122 + q2*fx2122 - fxy222 = p2*fx1222 + q2*fx2222 - - fxyz11 = p3*fxy111 + q3*fxy211 - fxyz21 = p3*fxy121 + q3*fxy221 - fxyz12 = p3*fxy112 + q3*fxy212 - fxyz22 = p3*fxy122 + q3*fxy222 - - fxyzq1 = p4*fxyz11 + q4*fxyz21 - fxyzq2 = p4*fxyz12 + q4*fxyz22 - - f = p5*fxyzq1 + q5*fxyzq2 + + f = p1*me%f(ix(1)) + q1*me%f(ix(2)) if (present(istat)) istat = 0 else if (present(istat)) istat = -1 + #:if NAME == 'dual' + allocate(f%der(size(x%der))) + #:endif f = zero end if - end subroutine interp_5d + end subroutine + #:endfor + + subroutine interp_1d_derivative(me, x, dfdx, istat) + use forwarddiff, only: derivative + class(linear_interp_1d),intent(inout) :: me + real(wp), intent(in) :: x + real(wp), intent(out) :: dfdx + integer, intent(out), optional :: istat + real(wp) :: f + call derivative(fcn, x, f, dfdx) + contains + function fcn(x_) result(res_) + use forwarddiff, only: dual + type(dual), intent(in) :: x_ + type(dual) :: res_ + call interp_1d_dual(me, x_, res_, istat) + end function + end subroutine + !***************************************************************************************** !***************************************************************************************** !> -! 6D linear interpolation routine. +! 2D linear interpolation routine. - pure subroutine interp_6d(me,x,y,z,q,r,s,f,istat) + pure subroutine interp_2d(me,x,y,f,istat) implicit none - class(linear_interp_6d),intent(inout) :: me + class(linear_interp_2d),intent(inout) :: me real(wp),intent(in) :: x real(wp),intent(in) :: y - real(wp),intent(in) :: z - real(wp),intent(in) :: q - real(wp),intent(in) :: r - real(wp),intent(in) :: s - real(wp),intent(out) :: f !! Interpolated \( f(x,y,z,q,r,s) \) - integer,intent(out),optional :: istat !! `0` : no problems, - !! `-1` : class has not been initialized - - integer,dimension(2) :: ix, iy, iz, iq, ir, is - real(wp) :: p1, p2, p3, p4, p5, p6 - real(wp) :: q1, q2, q3, q4, q5, q6 + real(wp),intent(out) :: f !! Interpolated \( f(x,y) \) + integer,intent(out),optional :: istat !! `0` : no problems, + !! `-1` : class has not been initialized + + integer,dimension(2) :: ix, iy + real(wp) :: p1, p2 + real(wp) :: q1, q2 integer :: mflag - real(wp) :: fx11111, fx21111, fx12111, fx22111, fx11211, fx21211, fx12211, & - fx22211, fxy1111, fxy2111, fxy1211, fxy2211, fxyz111, fxyz211, & - fxyzq11, fx11121, fx21121, fx12121, fx22121, fx11221, fx21221, & - fx12221, fx22221, fxy1121, fxy2121, fxy1221, fxy2221, fxyz121, & - fxyz221, fxyzq21, fx11112, fx21112, fx12112, fx22112, fx11212, & - fx21212, fx12212, fx22212, fxy1112, fxy2112, fxy1212, fxy2212, & - fxyz112, fxyz212, fxyzq12, fx11122, fx21122, fx12122, fx22122, & - fx11222, fx21222, fx12222, fx22222, fxy1122, fxy2122, fxy1222, & - fxy2222, fxyz122, fxyz222, fxyzq22, fxyzqr1, fxyzqr2 + real(wp) :: fx1, fx2 if (me%initialized) then call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag) call dintrv(me%y,y,me%iloy,iy(1),iy(2),mflag) - call dintrv(me%z,z,me%iloz,iz(1),iz(2),mflag) - call dintrv(me%q,q,me%iloq,iq(1),iq(2),mflag) - call dintrv(me%r,r,me%ilor,ir(1),ir(2),mflag) - call dintrv(me%s,s,me%ilos,is(1),is(2),mflag) q1 = (x-me%x(ix(1)))/(me%x(ix(2))-me%x(ix(1))) q2 = (y-me%y(iy(1)))/(me%y(iy(2))-me%y(iy(1))) - q3 = (z-me%z(iz(1)))/(me%z(iz(2))-me%z(iz(1))) - q4 = (q-me%q(iq(1)))/(me%q(iq(2))-me%q(iq(1))) - q5 = (r-me%r(ir(1)))/(me%r(ir(2))-me%r(ir(1))) - q6 = (s-me%s(is(1)))/(me%s(is(2))-me%s(is(1))) p1 = one-q1 p2 = one-q2 - p3 = one-q3 - p4 = one-q4 - p5 = one-q5 - p6 = one-q6 - - fx11111 = p1*me%f(ix(1),iy(1),iz(1),iq(1),ir(1),is(1)) + q1*me%f(ix(2),iy(1),iz(1),iq(1),ir(1),is(1)) - fx21111 = p1*me%f(ix(1),iy(2),iz(1),iq(1),ir(1),is(1)) + q1*me%f(ix(2),iy(2),iz(1),iq(1),ir(1),is(1)) - fx12111 = p1*me%f(ix(1),iy(1),iz(2),iq(1),ir(1),is(1)) + q1*me%f(ix(2),iy(1),iz(2),iq(1),ir(1),is(1)) - fx22111 = p1*me%f(ix(1),iy(2),iz(2),iq(1),ir(1),is(1)) + q1*me%f(ix(2),iy(2),iz(2),iq(1),ir(1),is(1)) - fx11211 = p1*me%f(ix(1),iy(1),iz(1),iq(2),ir(1),is(1)) + q1*me%f(ix(2),iy(1),iz(1),iq(2),ir(1),is(1)) - fx21211 = p1*me%f(ix(1),iy(2),iz(1),iq(2),ir(1),is(1)) + q1*me%f(ix(2),iy(2),iz(1),iq(2),ir(1),is(1)) - fx12211 = p1*me%f(ix(1),iy(1),iz(2),iq(2),ir(1),is(1)) + q1*me%f(ix(2),iy(1),iz(2),iq(2),ir(1),is(1)) - fx22211 = p1*me%f(ix(1),iy(2),iz(2),iq(2),ir(1),is(1)) + q1*me%f(ix(2),iy(2),iz(2),iq(2),ir(1),is(1)) - fx11121 = p1*me%f(ix(1),iy(1),iz(1),iq(1),ir(2),is(1)) + q1*me%f(ix(2),iy(1),iz(1),iq(1),ir(2),is(1)) - fx21121 = p1*me%f(ix(1),iy(2),iz(1),iq(1),ir(2),is(1)) + q1*me%f(ix(2),iy(2),iz(1),iq(1),ir(2),is(1)) - fx12121 = p1*me%f(ix(1),iy(1),iz(2),iq(1),ir(2),is(1)) + q1*me%f(ix(2),iy(1),iz(2),iq(1),ir(2),is(1)) - fx22121 = p1*me%f(ix(1),iy(2),iz(2),iq(1),ir(2),is(1)) + q1*me%f(ix(2),iy(2),iz(2),iq(1),ir(2),is(1)) - fx11221 = p1*me%f(ix(1),iy(1),iz(1),iq(2),ir(2),is(1)) + q1*me%f(ix(2),iy(1),iz(1),iq(2),ir(2),is(1)) - fx21221 = p1*me%f(ix(1),iy(2),iz(1),iq(2),ir(2),is(1)) + q1*me%f(ix(2),iy(2),iz(1),iq(2),ir(2),is(1)) - fx12221 = p1*me%f(ix(1),iy(1),iz(2),iq(2),ir(2),is(1)) + q1*me%f(ix(2),iy(1),iz(2),iq(2),ir(2),is(1)) - fx22221 = p1*me%f(ix(1),iy(2),iz(2),iq(2),ir(2),is(1)) + q1*me%f(ix(2),iy(2),iz(2),iq(2),ir(2),is(1)) - fx11112 = p1*me%f(ix(1),iy(1),iz(1),iq(1),ir(1),is(2)) + q1*me%f(ix(2),iy(1),iz(1),iq(1),ir(1),is(2)) - fx21112 = p1*me%f(ix(1),iy(2),iz(1),iq(1),ir(1),is(2)) + q1*me%f(ix(2),iy(2),iz(1),iq(1),ir(1),is(2)) - fx12112 = p1*me%f(ix(1),iy(1),iz(2),iq(1),ir(1),is(2)) + q1*me%f(ix(2),iy(1),iz(2),iq(1),ir(1),is(2)) - fx22112 = p1*me%f(ix(1),iy(2),iz(2),iq(1),ir(1),is(2)) + q1*me%f(ix(2),iy(2),iz(2),iq(1),ir(1),is(2)) - fx11212 = p1*me%f(ix(1),iy(1),iz(1),iq(2),ir(1),is(2)) + q1*me%f(ix(2),iy(1),iz(1),iq(2),ir(1),is(2)) - fx21212 = p1*me%f(ix(1),iy(2),iz(1),iq(2),ir(1),is(2)) + q1*me%f(ix(2),iy(2),iz(1),iq(2),ir(1),is(2)) - fx12212 = p1*me%f(ix(1),iy(1),iz(2),iq(2),ir(1),is(2)) + q1*me%f(ix(2),iy(1),iz(2),iq(2),ir(1),is(2)) - fx22212 = p1*me%f(ix(1),iy(2),iz(2),iq(2),ir(1),is(2)) + q1*me%f(ix(2),iy(2),iz(2),iq(2),ir(1),is(2)) - fx11122 = p1*me%f(ix(1),iy(1),iz(1),iq(1),ir(2),is(2)) + q1*me%f(ix(2),iy(1),iz(1),iq(1),ir(2),is(2)) - fx21122 = p1*me%f(ix(1),iy(2),iz(1),iq(1),ir(2),is(2)) + q1*me%f(ix(2),iy(2),iz(1),iq(1),ir(2),is(2)) - fx12122 = p1*me%f(ix(1),iy(1),iz(2),iq(1),ir(2),is(2)) + q1*me%f(ix(2),iy(1),iz(2),iq(1),ir(2),is(2)) - fx22122 = p1*me%f(ix(1),iy(2),iz(2),iq(1),ir(2),is(2)) + q1*me%f(ix(2),iy(2),iz(2),iq(1),ir(2),is(2)) - fx11222 = p1*me%f(ix(1),iy(1),iz(1),iq(2),ir(2),is(2)) + q1*me%f(ix(2),iy(1),iz(1),iq(2),ir(2),is(2)) - fx21222 = p1*me%f(ix(1),iy(2),iz(1),iq(2),ir(2),is(2)) + q1*me%f(ix(2),iy(2),iz(1),iq(2),ir(2),is(2)) - fx12222 = p1*me%f(ix(1),iy(1),iz(2),iq(2),ir(2),is(2)) + q1*me%f(ix(2),iy(1),iz(2),iq(2),ir(2),is(2)) - fx22222 = p1*me%f(ix(1),iy(2),iz(2),iq(2),ir(2),is(2)) + q1*me%f(ix(2),iy(2),iz(2),iq(2),ir(2),is(2)) - - fxy1111 = p2*fx11111 + q2*fx21111 - fxy2111 = p2*fx12111 + q2*fx22111 - fxy1211 = p2*fx11211 + q2*fx21211 - fxy2211 = p2*fx12211 + q2*fx22211 - fxy1121 = p2*fx11121 + q2*fx21121 - fxy2121 = p2*fx12121 + q2*fx22121 - fxy1221 = p2*fx11221 + q2*fx21221 - fxy2221 = p2*fx12221 + q2*fx22221 - fxy1112 = p2*fx11112 + q2*fx21112 - fxy2112 = p2*fx12112 + q2*fx22112 - fxy1212 = p2*fx11212 + q2*fx21212 - fxy2212 = p2*fx12212 + q2*fx22212 - fxy1122 = p2*fx11122 + q2*fx21122 - fxy2122 = p2*fx12122 + q2*fx22122 - fxy1222 = p2*fx11222 + q2*fx21222 - fxy2222 = p2*fx12222 + q2*fx22222 - - fxyz111 = p3*fxy1111 + q3*fxy2111 - fxyz211 = p3*fxy1211 + q3*fxy2211 - fxyz121 = p3*fxy1121 + q3*fxy2121 - fxyz221 = p3*fxy1221 + q3*fxy2221 - fxyz112 = p3*fxy1112 + q3*fxy2112 - fxyz212 = p3*fxy1212 + q3*fxy2212 - fxyz122 = p3*fxy1122 + q3*fxy2122 - fxyz222 = p3*fxy1222 + q3*fxy2222 - - fxyzq11 = p4*fxyz111 + q4*fxyz211 - fxyzq21 = p4*fxyz121 + q4*fxyz221 - fxyzq12 = p4*fxyz112 + q4*fxyz212 - fxyzq22 = p4*fxyz122 + q4*fxyz222 - - fxyzqr1 = p5*fxyzq11 + q5*fxyzq21 - fxyzqr2 = p5*fxyzq12 + q5*fxyzq22 - - f = p6*fxyzqr1 + q6*fxyzqr2 + + fx1 = p1*me%f(ix(1),iy(1)) + q1*me%f(ix(2),iy(1)) + fx2 = p1*me%f(ix(1),iy(2)) + q1*me%f(ix(2),iy(2)) + + f = p2*fx1 + q2*fx2 if (present(istat)) istat = 0 else @@ -1154,7 +334,7 @@ pure subroutine interp_6d(me,x,y,z,q,r,s,f,istat) end if - end subroutine interp_6d + end subroutine interp_2d !***************************************************************************************** !***************************************************************************************** @@ -1180,13 +360,16 @@ end subroutine interp_6d ! * Jacob Williams, 2/22/2016 : modified bspline-fortran `dintrv` routine for ! linear interpolation/extrapolation use. ! * Jacob Williams, 10/9/2019 : added optional `inearest` output. - - pure subroutine dintrv(xt,x,ilo,ileft,iright,mflag,inearest) + #:for TYPE1, NAME in TYPES_NAMES + pure subroutine dintrv_${NAME}$(xt,x,ilo,ileft,iright,mflag,inearest) + #:if NAME == 'dual' + use forwarddiff + #:endif implicit none real(wp),dimension(:),intent(in) :: xt !! a knot or break point vector - real(wp),intent(in) :: x !! argument + ${TYPE1}$,intent(in) :: x !! argument integer,intent(inout) :: ilo !! an initialization parameter which must be set !! to 1 the first time the array `xt` is !! processed by dintrv. `ilo` contains information for @@ -1315,247 +498,8 @@ pure subroutine dintrv(xt,x,ilo,ileft,iright,mflag,inearest) endif end do - end subroutine dintrv -!***************************************************************************************** - -!***************************************************************************************** -!> -! 1D nearest neighbor interpolation routine. - - pure subroutine nearest_1d(me,x,f,istat) - - implicit none - - class(nearest_interp_1d),intent(inout) :: me - real(wp),intent(in) :: x - real(wp),intent(out) :: f !! Nearest \( f(x) \) - integer,intent(out),optional :: istat !! `0` : no problems, - !! `-1` : class has not been initialized - - integer :: mflag - integer,dimension(2) :: ix - integer :: i - - if (me%initialized) then - - call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag,i) - - f = me%f(i) - if (present(istat)) istat = 0 - - else - - if (present(istat)) istat = -1 - f = zero - - end if - - end subroutine nearest_1d -!***************************************************************************************** - -!***************************************************************************************** -!> -! 2D nearest neighbor interpolation routine. - - pure subroutine nearest_2d(me,x,y,f,istat) - - implicit none - - class(nearest_interp_2d),intent(inout) :: me - real(wp),intent(in) :: x - real(wp),intent(in) :: y - real(wp),intent(out) :: f !! Nearest \( f(x,y) \) - integer,intent(out),optional :: istat !! `0` : no problems, - !! `-1` : class has not been initialized - - integer :: mflag - integer,dimension(2) :: ix, iy - integer :: i, j - - if (me%initialized) then - - call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag,i) - call dintrv(me%y,y,me%iloy,iy(1),iy(2),mflag,j) - - f = me%f(i,j) - if (present(istat)) istat = 0 - - else - - if (present(istat)) istat = -1 - f = zero - - end if - - end subroutine nearest_2d -!***************************************************************************************** - -!***************************************************************************************** -!> -! 3D nearest neighbor interpolation routine. - - pure subroutine nearest_3d(me,x,y,z,f,istat) - - implicit none - - class(nearest_interp_3d),intent(inout) :: me - real(wp),intent(in) :: x - real(wp),intent(in) :: y - real(wp),intent(in) :: z - real(wp),intent(out) :: f !! Nearest \( f(x,y,z) \) - integer,intent(out),optional :: istat !! `0` : no problems, - !! `-1` : class has not been initialized - - integer :: mflag - integer,dimension(2) :: ix, iy, iz - integer :: i, j, k - - if (me%initialized) then - - call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag,i) - call dintrv(me%y,y,me%iloy,iy(1),iy(2),mflag,j) - call dintrv(me%z,z,me%iloz,iz(1),iz(2),mflag,k) - - f = me%f(i,j,k) - if (present(istat)) istat = 0 - - else - - if (present(istat)) istat = -1 - f = zero - - end if - - end subroutine nearest_3d -!***************************************************************************************** - -!***************************************************************************************** -!> -! 4D nearest neighbor interpolation routine. - - pure subroutine nearest_4d(me,x,y,z,q,f,istat) - - implicit none - - class(nearest_interp_4d),intent(inout) :: me - real(wp),intent(in) :: x - real(wp),intent(in) :: y - real(wp),intent(in) :: z - real(wp),intent(in) :: q - real(wp),intent(out) :: f !! Nearest \( f(x,y,z,q) \) - integer,intent(out),optional :: istat !! `0` : no problems, - !! `-1` : class has not been initialized - - integer :: mflag - integer,dimension(2) :: ix, iy, iz, iq - integer :: i, j, k, l - - if (me%initialized) then - - call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag,i) - call dintrv(me%y,y,me%iloy,iy(1),iy(2),mflag,j) - call dintrv(me%z,z,me%iloz,iz(1),iz(2),mflag,k) - call dintrv(me%q,q,me%iloq,iq(1),iq(2),mflag,l) - - f = me%f(i,j,k,l) - if (present(istat)) istat = 0 - - else - - if (present(istat)) istat = -1 - f = zero - - end if - - end subroutine nearest_4d -!***************************************************************************************** - -!***************************************************************************************** -!> -! 5D nearest neighbor interpolation routine. - - pure subroutine nearest_5d(me,x,y,z,q,r,f,istat) - - implicit none - - class(nearest_interp_5d),intent(inout) :: me - real(wp),intent(in) :: x - real(wp),intent(in) :: y - real(wp),intent(in) :: z - real(wp),intent(in) :: q - real(wp),intent(in) :: r - real(wp),intent(out) :: f !! Nearest \( f(x,y,z,q,r) \) - integer,intent(out),optional :: istat !! `0` : no problems, - !! `-1` : class has not been initialized - - integer :: mflag - integer,dimension(2) :: ix, iy, iz, iq, ir - integer :: i, j, k, l, m - - if (me%initialized) then - - call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag,i) - call dintrv(me%y,y,me%iloy,iy(1),iy(2),mflag,j) - call dintrv(me%z,z,me%iloz,iz(1),iz(2),mflag,k) - call dintrv(me%q,q,me%iloq,iq(1),iq(2),mflag,l) - call dintrv(me%r,r,me%ilor,ir(1),ir(2),mflag,m) - - f = me%f(i,j,k,l,m) - if (present(istat)) istat = 0 - - else - - if (present(istat)) istat = -1 - f = zero - - end if - - end subroutine nearest_5d -!***************************************************************************************** - -!***************************************************************************************** -!> -! 6D nearest neighbor interpolation routine. - - pure subroutine nearest_6d(me,x,y,z,q,r,s,f,istat) - - implicit none - - class(nearest_interp_6d),intent(inout) :: me - real(wp),intent(in) :: x - real(wp),intent(in) :: y - real(wp),intent(in) :: z - real(wp),intent(in) :: q - real(wp),intent(in) :: r - real(wp),intent(in) :: s - real(wp),intent(out) :: f !! Nearest \( f(x,y,z,q,r,s) \) - integer,intent(out),optional :: istat !! `0` : no problems, - !! `-1` : class has not been initialized - - integer :: mflag - integer,dimension(2) :: ix, iy, iz, iq, ir, is - integer :: i, j, k, l, m, n - - if (me%initialized) then - - call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag,i) - call dintrv(me%y,y,me%iloy,iy(1),iy(2),mflag,j) - call dintrv(me%z,z,me%iloz,iz(1),iz(2),mflag,k) - call dintrv(me%q,q,me%iloq,iq(1),iq(2),mflag,l) - call dintrv(me%r,r,me%ilor,ir(1),ir(2),mflag,m) - call dintrv(me%s,s,me%ilos,is(1),is(2),mflag,n) - - f = me%f(i,j,k,l,m,n) - if (present(istat)) istat = 0 - - else - - if (present(istat)) istat = -1 - f = zero - - end if - - end subroutine nearest_6d + end subroutine + #:endfor !***************************************************************************************** !***************************************************************************************** @@ -1601,12 +545,6 @@ pure subroutine check_inputs(me,x,y,z,q,r,s,ierr) if (ierr == 0) then select type (me) - class is (nearest_interp_1d) - class is (nearest_interp_2d) - class is (nearest_interp_3d) - class is (nearest_interp_4d) - class is (nearest_interp_5d) - class is (nearest_interp_6d) class default ! need at least two points for linear interpolation: if (size(x)==1) ierr = 100