!***************************************************************************************** !> author: Jacob Williams ! license: BSD ! date: 12/6/2015 ! ! Object-oriented style wrappers to [[bspline_sub_module]]. ! This module provides classes ([[bspline_1d(type)]], [[bspline_2d(type)]], ! [[bspline_3d(type)]], [[bspline_4d(type)]], [[bspline_5d(type)]], and [[bspline_6d(type)]]) ! which can be used instead of the main subroutine interface. module bspline_oo_module use bspline_kinds_module, only: wp, ip use,intrinsic :: iso_fortran_env, only: error_unit use bspline_sub_module implicit none private integer(ip),parameter :: int_size = storage_size(1_ip,kind=ip) !! size of a default integer [bits] integer(ip),parameter :: logical_size = storage_size(.true.,kind=ip) !! size of a default logical [bits] integer(ip),parameter :: real_size = storage_size(1.0_wp,kind=ip) !! size of a `real(wp)` [bits] type,public,abstract :: bspline_class !! Base class for the b-spline types private integer(ip) :: inbvx = 1_ip !! internal variable used by [[dbvalu]] for efficient processing integer(ip) :: iflag = 1_ip !! saved `iflag` from the list routine call. logical :: initialized = .false. !! true if the class is initialized and ready to use logical :: extrap = .false. !! if true, then extrapolation is allowed during evaluation integer(ip) :: nf = 0_ip !! number of outputs (the last dimension of `bcoef`) contains private procedure,non_overridable :: destroy_base !! destructor for the abstract type procedure,non_overridable :: set_extrap_flag !! internal routine to set the `extrap` flag procedure(destroy_func),deferred,public :: destroy !! destructor procedure(size_func),deferred,public :: size_of !! size of the structure in bits procedure,public,non_overridable :: status_ok !! returns true if the last `iflag` status code was `=0`. procedure,public,non_overridable :: status_message => get_bspline_status_message !! retrieve the last !! status message procedure,public,non_overridable :: clear_flag => clear_bspline_flag !! to reset the `iflag` saved in the class. procedure,public,non_overridable :: number_of_outputs !! number of outputs (functions) of the spline end type bspline_class abstract interface pure subroutine destroy_func(me) !! interface for bspline destructor routines import :: bspline_class implicit none class(bspline_class),intent(inout) :: me end subroutine destroy_func pure function size_func(me) result(s) !! interface for size routines import :: bspline_class,ip implicit none class(bspline_class),intent(in) :: me integer(ip) :: s !! size of the structure in bits end function size_func end interface type,extends(bspline_class),public :: bspline_1d !! Class for 1d b-spline interpolation. !! !!@note The 1D class also contains methods for computing !! definite integrals, roots, the minimum and maximum, !! and the control points. private integer(ip) :: nx = 0_ip !! Number of \(x\) abcissae integer(ip) :: kx = 0_ip !! The order of spline pieces in \(x\) real(wp),dimension(:,:),allocatable :: bcoef !! array of coefficients of the b-spline interpolant !! (the last dimension is the output index) real(wp),dimension(:),allocatable :: tx !! The knots in the \(x\) direction for the spline interpolant real(wp),dimension(:),allocatable :: work_val_1 !! [[db1val] work array of dimension `3*kx` contains private generic,public :: initialize => initialize_1d_auto_knots,initialize_1d_specify_knots,& initialize_1d_auto_knots_multi,initialize_1d_specify_knots_multi procedure :: initialize_1d_auto_knots_multi procedure :: initialize_1d_specify_knots_multi procedure :: initialize_1d_auto_knots procedure :: initialize_1d_specify_knots generic,public :: evaluate => evaluate_1d, evaluate_1d_multi procedure :: evaluate_1d procedure :: evaluate_1d_multi procedure,public :: destroy => destroy_1d procedure,public :: size_of => size_1d generic,public :: initialize_from_coefficients => initialize_from_coefficients_1d,& initialize_from_coefficients_1d_multi procedure :: initialize_from_coefficients_1d procedure :: initialize_from_coefficients_1d_multi procedure,public :: get_knots => get_knots_1d generic,public :: get_coefficients => get_coefficients_1d, get_coefficients_1d_multi procedure :: get_coefficients_1d procedure :: get_coefficients_1d_multi procedure,public :: control_points => control_points_1d procedure,public :: integral => integral_1d procedure,public :: fintegral => fintegral_1d procedure,public :: roots => roots_1d procedure,public :: inverse => inverse_1d procedure,public :: minimum => minimum_1d procedure,public :: maximum => maximum_1d final :: finalize_1d end type bspline_1d type,extends(bspline_class),public :: bspline_2d !! Class for 2d b-spline interpolation. private integer(ip) :: nx = 0_ip !! Number of \(x\) abcissae integer(ip) :: ny = 0_ip !! Number of \(y\) abcissae integer(ip) :: kx = 0_ip !! The order of spline pieces in \(x\) integer(ip) :: ky = 0_ip !! The order of spline pieces in \(y\) real(wp),dimension(:,:,:),allocatable :: bcoef !! array of coefficients of the b-spline interpolant !! (the last dimension is the output index) real(wp),dimension(:),allocatable :: tx !! The knots in the \(x\) direction for the spline interpolant real(wp),dimension(:),allocatable :: ty !! The knots in the \(y\) direction for the spline interpolant integer(ip) :: inbvy = 1_ip !! internal variable used for efficient processing integer(ip) :: iloy = 1_ip !! internal variable used for efficient processing real(wp),dimension(:),allocatable :: work_val_1 !! [[db2val] work array of dimension `ky` real(wp),dimension(:),allocatable :: work_val_2 !! [[db2val] work array of dimension `3_ip*max(kx,ky)` contains private generic,public :: initialize => initialize_2d_auto_knots,initialize_2d_specify_knots,& initialize_2d_auto_knots_multi,initialize_2d_specify_knots_multi procedure :: initialize_2d_auto_knots_multi procedure :: initialize_2d_specify_knots_multi procedure :: initialize_2d_auto_knots procedure :: initialize_2d_specify_knots generic,public :: evaluate => evaluate_2d, evaluate_2d_multi procedure :: evaluate_2d procedure :: evaluate_2d_multi procedure,public :: destroy => destroy_2d procedure,public :: size_of => size_2d generic,public :: initialize_from_coefficients => initialize_from_coefficients_2d,& initialize_from_coefficients_2d_multi procedure :: initialize_from_coefficients_2d procedure :: initialize_from_coefficients_2d_multi procedure,public :: get_knots => get_knots_2d generic,public :: get_coefficients => get_coefficients_2d, get_coefficients_2d_multi procedure :: get_coefficients_2d procedure :: get_coefficients_2d_multi final :: finalize_2d end type bspline_2d type,extends(bspline_class),public :: bspline_3d !! Class for 3d b-spline interpolation. private integer(ip) :: nx = 0_ip !! Number of \(x\) abcissae integer(ip) :: ny = 0_ip !! Number of \(y\) abcissae integer(ip) :: nz = 0_ip !! Number of \(z\) abcissae integer(ip) :: kx = 0_ip !! The order of spline pieces in \(x\) integer(ip) :: ky = 0_ip !! The order of spline pieces in \(y\) integer(ip) :: kz = 0_ip !! The order of spline pieces in \(z\) real(wp),dimension(:,:,:,:),allocatable :: bcoef !! array of coefficients of the b-spline interpolant !! (the last dimension is the output index) real(wp),dimension(:),allocatable :: tx !! The knots in the \(x\) direction for the spline interpolant real(wp),dimension(:),allocatable :: ty !! The knots in the \(y\) direction for the spline interpolant real(wp),dimension(:),allocatable :: tz !! The knots in the \(z\) direction for the spline interpolant integer(ip) :: inbvy = 1_ip !! internal variable used for efficient processing integer(ip) :: inbvz = 1_ip !! internal variable used for efficient processing integer(ip) :: iloy = 1_ip !! internal variable used for efficient processing integer(ip) :: iloz = 1_ip !! internal variable used for efficient processing real(wp),dimension(:,:),allocatable :: work_val_1 !! [[db3val] work array of dimension `ky,kz` real(wp),dimension(:),allocatable :: work_val_2 !! [[db3val] work array of dimension `kz` real(wp),dimension(:),allocatable :: work_val_3 !! [[db3val] work array of dimension `3_ip*max(kx,ky,kz)` contains private generic,public :: initialize => initialize_3d_auto_knots,initialize_3d_specify_knots,& initialize_3d_auto_knots_multi,initialize_3d_specify_knots_multi procedure :: initialize_3d_auto_knots_multi procedure :: initialize_3d_specify_knots_multi procedure :: initialize_3d_auto_knots procedure :: initialize_3d_specify_knots generic,public :: evaluate => evaluate_3d, evaluate_3d_multi procedure :: evaluate_3d procedure :: evaluate_3d_multi procedure,public :: destroy => destroy_3d procedure,public :: size_of => size_3d generic,public :: initialize_from_coefficients => initialize_from_coefficients_3d,& initialize_from_coefficients_3d_multi procedure :: initialize_from_coefficients_3d procedure :: initialize_from_coefficients_3d_multi procedure,public :: get_knots => get_knots_3d generic,public :: get_coefficients => get_coefficients_3d, get_coefficients_3d_multi procedure :: get_coefficients_3d procedure :: get_coefficients_3d_multi final :: finalize_3d end type bspline_3d type,extends(bspline_class),public :: bspline_4d !! Class for 4d b-spline interpolation. private integer(ip) :: nx = 0_ip !! Number of \(x\) abcissae integer(ip) :: ny = 0_ip !! Number of \(y\) abcissae integer(ip) :: nz = 0_ip !! Number of \(z\) abcissae integer(ip) :: nq = 0_ip !! Number of \(q\) abcissae integer(ip) :: kx = 0_ip !! The order of spline pieces in \(x\) integer(ip) :: ky = 0_ip !! The order of spline pieces in \(y\) integer(ip) :: kz = 0_ip !! The order of spline pieces in \(z\) integer(ip) :: kq = 0_ip !! The order of spline pieces in \(q\) real(wp),dimension(:,:,:,:,:),allocatable :: bcoef !! array of coefficients of the b-spline interpolant !! (the last dimension is the output index) real(wp),dimension(:),allocatable :: tx !! The knots in the \(x\) direction for the spline interpolant real(wp),dimension(:),allocatable :: ty !! The knots in the \(y\) direction for the spline interpolant real(wp),dimension(:),allocatable :: tz !! The knots in the \(z\) direction for the spline interpolant real(wp),dimension(:),allocatable :: tq !! The knots in the \(q\) direction for the spline interpolant integer(ip) :: inbvy = 1_ip !! internal variable used for efficient processing integer(ip) :: inbvz = 1_ip !! internal variable used for efficient processing integer(ip) :: inbvq = 1_ip !! internal variable used for efficient processing integer(ip) :: iloy = 1_ip !! internal variable used for efficient processing integer(ip) :: iloz = 1_ip !! internal variable used for efficient processing integer(ip) :: iloq = 1_ip !! internal variable used for efficient processing real(wp),dimension(:,:,:),allocatable :: work_val_1 !! [[db4val]] work array of dimension `ky,kz,kq` real(wp),dimension(:,:),allocatable :: work_val_2 !! [[db4val]] work array of dimension `kz,kq` real(wp),dimension(:),allocatable :: work_val_3 !! [[db4val]] work array of dimension `kq` real(wp),dimension(:),allocatable :: work_val_4 !! [[db4val]] work array of dimension `3_ip*max(kx,ky,kz,kq)` contains private generic,public :: initialize => initialize_4d_auto_knots,initialize_4d_specify_knots,& initialize_4d_auto_knots_multi,initialize_4d_specify_knots_multi procedure :: initialize_4d_auto_knots_multi procedure :: initialize_4d_specify_knots_multi procedure :: initialize_4d_auto_knots procedure :: initialize_4d_specify_knots generic,public :: evaluate => evaluate_4d, evaluate_4d_multi procedure :: evaluate_4d procedure :: evaluate_4d_multi procedure,public :: destroy => destroy_4d procedure,public :: size_of => size_4d generic,public :: initialize_from_coefficients => initialize_from_coefficients_4d,& initialize_from_coefficients_4d_multi procedure :: initialize_from_coefficients_4d procedure :: initialize_from_coefficients_4d_multi procedure,public :: get_knots => get_knots_4d generic,public :: get_coefficients => get_coefficients_4d, get_coefficients_4d_multi procedure :: get_coefficients_4d procedure :: get_coefficients_4d_multi final :: finalize_4d end type bspline_4d type,extends(bspline_class),public :: bspline_5d !! Class for 5d b-spline interpolation. private integer(ip) :: nx = 0_ip !! Number of \(x\) abcissae integer(ip) :: ny = 0_ip !! Number of \(y\) abcissae integer(ip) :: nz = 0_ip !! Number of \(z\) abcissae integer(ip) :: nq = 0_ip !! Number of \(q\) abcissae integer(ip) :: nr = 0_ip !! Number of \(r\) abcissae integer(ip) :: kx = 0_ip !! The order of spline pieces in \(x\) integer(ip) :: ky = 0_ip !! The order of spline pieces in \(y\) integer(ip) :: kz = 0_ip !! The order of spline pieces in \(z\) integer(ip) :: kq = 0_ip !! The order of spline pieces in \(q\) integer(ip) :: kr = 0_ip !! The order of spline pieces in \(r\) real(wp),dimension(:,:,:,:,:,:),allocatable :: bcoef !! array of coefficients of the b-spline interpolant !! (the last dimension is the output index) real(wp),dimension(:),allocatable :: tx !! The knots in the \(x\) direction for the spline interpolant real(wp),dimension(:),allocatable :: ty !! The knots in the \(y\) direction for the spline interpolant real(wp),dimension(:),allocatable :: tz !! The knots in the \(z\) direction for the spline interpolant real(wp),dimension(:),allocatable :: tq !! The knots in the \(q\) direction for the spline interpolant real(wp),dimension(:),allocatable :: tr !! The knots in the \(r\) direction for the spline interpolant integer(ip) :: inbvy = 1_ip !! internal variable used for efficient processing integer(ip) :: inbvz = 1_ip !! internal variable used for efficient processing integer(ip) :: inbvq = 1_ip !! internal variable used for efficient processing integer(ip) :: inbvr = 1_ip !! internal variable used for efficient processing integer(ip) :: iloy = 1_ip !! internal variable used for efficient processing integer(ip) :: iloz = 1_ip !! internal variable used for efficient processing integer(ip) :: iloq = 1_ip !! internal variable used for efficient processing integer(ip) :: ilor = 1_ip !! internal variable used for efficient processing real(wp),dimension(:,:,:,:),allocatable :: work_val_1 !! [[db5val]] work array of dimension `ky,kz,kq,kr` real(wp),dimension(:,:,:),allocatable :: work_val_2 !! [[db5val]] work array of dimension `kz,kq,kr` real(wp),dimension(:,:),allocatable :: work_val_3 !! [[db5val]] work array of dimension `kq,kr` real(wp),dimension(:),allocatable :: work_val_4 !! [[db5val]] work array of dimension `kr` real(wp),dimension(:),allocatable :: work_val_5 !! [[db5val]] work array of dimension `3_ip*max(kx,ky,kz,kq,kr)` contains private generic,public :: initialize => initialize_5d_auto_knots,initialize_5d_specify_knots,& initialize_5d_auto_knots_multi,initialize_5d_specify_knots_multi procedure :: initialize_5d_auto_knots_multi procedure :: initialize_5d_specify_knots_multi procedure :: initialize_5d_auto_knots procedure :: initialize_5d_specify_knots generic,public :: evaluate => evaluate_5d, evaluate_5d_multi procedure :: evaluate_5d procedure :: evaluate_5d_multi procedure,public :: destroy => destroy_5d procedure,public :: size_of => size_5d generic,public :: initialize_from_coefficients => initialize_from_coefficients_5d,& initialize_from_coefficients_5d_multi procedure :: initialize_from_coefficients_5d procedure :: initialize_from_coefficients_5d_multi procedure,public :: get_knots => get_knots_5d generic,public :: get_coefficients => get_coefficients_5d, get_coefficients_5d_multi procedure :: get_coefficients_5d procedure :: get_coefficients_5d_multi final :: finalize_5d end type bspline_5d type,extends(bspline_class),public :: bspline_6d !! Class for 6d b-spline interpolation. private integer(ip) :: nx = 0_ip !! Number of \(x\) abcissae integer(ip) :: ny = 0_ip !! Number of \(y\) abcissae integer(ip) :: nz = 0_ip !! Number of \(z\) abcissae integer(ip) :: nq = 0_ip !! Number of \(q\) abcissae integer(ip) :: nr = 0_ip !! Number of \(r\) abcissae integer(ip) :: ns = 0_ip !! Number of \(s\) abcissae integer(ip) :: kx = 0_ip !! The order of spline pieces in \(x\) integer(ip) :: ky = 0_ip !! The order of spline pieces in \(y\) integer(ip) :: kz = 0_ip !! The order of spline pieces in \(z\) integer(ip) :: kq = 0_ip !! The order of spline pieces in \(q\) integer(ip) :: kr = 0_ip !! The order of spline pieces in \(r\) integer(ip) :: ks = 0_ip !! The order of spline pieces in \(s\) real(wp),dimension(:,:,:,:,:,:,:),allocatable :: bcoef !! array of coefficients of the b-spline interpolant !! (the last dimension is the output index) real(wp),dimension(:),allocatable :: tx !! The knots in the \(x\) direction for the spline interpolant real(wp),dimension(:),allocatable :: ty !! The knots in the \(y\) direction for the spline interpolant real(wp),dimension(:),allocatable :: tz !! The knots in the \(z\) direction for the spline interpolant real(wp),dimension(:),allocatable :: tq !! The knots in the \(q\) direction for the spline interpolant real(wp),dimension(:),allocatable :: tr !! The knots in the \(r\) direction for the spline interpolant real(wp),dimension(:),allocatable :: ts !! The knots in the \(s\) direction for the spline interpolant integer(ip) :: inbvy = 1_ip !! internal variable used for efficient processing integer(ip) :: inbvz = 1_ip !! internal variable used for efficient processing integer(ip) :: inbvq = 1_ip !! internal variable used for efficient processing integer(ip) :: inbvr = 1_ip !! internal variable used for efficient processing integer(ip) :: inbvs = 1_ip !! internal variable used for efficient processing integer(ip) :: iloy = 1_ip !! internal variable used for efficient processing integer(ip) :: iloz = 1_ip !! internal variable used for efficient processing integer(ip) :: iloq = 1_ip !! internal variable used for efficient processing integer(ip) :: ilor = 1_ip !! internal variable used for efficient processing integer(ip) :: ilos = 1_ip !! internal variable used for efficient processing real(wp),dimension(:,:,:,:,:),allocatable :: work_val_1 !! [[db6val]] work array of dimension `ky,kz,kq,kr,ks` real(wp),dimension(:,:,:,:),allocatable :: work_val_2 !! [[db6val]] work array of dimension `kz,kq,kr,ks` real(wp),dimension(:,:,:),allocatable :: work_val_3 !! [[db6val]] work array of dimension `kq,kr,ks` real(wp),dimension(:,:),allocatable :: work_val_4 !! [[db6val]] work array of dimension `kr,ks` real(wp),dimension(:),allocatable :: work_val_5 !! [[db6val]] work array of dimension `ks` real(wp),dimension(:),allocatable :: work_val_6 !! [[db6val]] work array of dimension `3_ip*max(kx,ky,kz,kq,kr,ks)` contains private generic,public :: initialize => initialize_6d_auto_knots,initialize_6d_specify_knots,& initialize_6d_auto_knots_multi,initialize_6d_specify_knots_multi procedure :: initialize_6d_auto_knots_multi procedure :: initialize_6d_specify_knots_multi procedure :: initialize_6d_auto_knots procedure :: initialize_6d_specify_knots generic,public :: evaluate => evaluate_6d, evaluate_6d_multi procedure :: evaluate_6d procedure :: evaluate_6d_multi procedure,public :: destroy => destroy_6d procedure,public :: size_of => size_6d generic,public :: initialize_from_coefficients => initialize_from_coefficients_6d,& initialize_from_coefficients_6d_multi procedure :: initialize_from_coefficients_6d procedure :: initialize_from_coefficients_6d_multi procedure,public :: get_knots => get_knots_6d generic,public :: get_coefficients => get_coefficients_6d, get_coefficients_6d_multi procedure :: get_coefficients_6d procedure :: get_coefficients_6d_multi final :: finalize_6d end type bspline_6d interface bspline_1d !! Constructor for [[bspline_1d(type)]] procedure :: bspline_1d_constructor_empty,& bspline_1d_constructor_auto_knots,& bspline_1d_constructor_specify_knots end interface interface bspline_2d !! Constructor for [[bspline_2d(type)]] procedure :: bspline_2d_constructor_empty,& bspline_2d_constructor_auto_knots,& bspline_2d_constructor_specify_knots end interface interface bspline_3d !! Constructor for [[bspline_3d(type)]] procedure :: bspline_3d_constructor_empty,& bspline_3d_constructor_auto_knots,& bspline_3d_constructor_specify_knots end interface interface bspline_4d !! Constructor for [[bspline_4d(type)]] procedure :: bspline_4d_constructor_empty,& bspline_4d_constructor_auto_knots,& bspline_4d_constructor_specify_knots end interface interface bspline_5d !! Constructor for [[bspline_5d(type)]] procedure :: bspline_5d_constructor_empty,& bspline_5d_constructor_auto_knots,& bspline_5d_constructor_specify_knots end interface interface bspline_6d !! Constructor for [[bspline_6d(type)]] procedure :: bspline_6d_constructor_empty,& bspline_6d_constructor_auto_knots,& bspline_6d_constructor_specify_knots end interface contains !***************************************************************************************** !***************************************************************************************** !> ! Returns the number of outputs (functions of the same inputs) of the ! spline: 1 for a spline initialized with a single function, and `nf` ! for one initialized with `nf` functions (the last dimension of `fcn`). ! Returns 0 if the class is not initialized. pure function number_of_outputs(me) result(nf) implicit none class(bspline_class),intent(in) :: me !! the spline integer(ip) :: nf !! number of outputs nf = me%nf end function number_of_outputs !***************************************************************************************** !***************************************************************************************** !> ! This routines returns true if the `iflag` code from the last ! routine called was `=0`. Maybe of the routines have output `iflag` ! variables, so they can be checked explicitly, or this routine ! can be used. ! ! If the class is initialized using a function constructor, then ! this is the only way to know if it was properly initialized, ! since those are pure functions with not output `iflag` arguments. ! ! If `status_ok=.false.`, then the error message can be ! obtained from the [[get_bspline_status_message]] routine. ! ! Note: after an error condition, the [[clear_bspline_flag]] routine ! can be called to reset the `iflag` to 0. elemental function status_ok(me) result(ok) implicit none class(bspline_class),intent(in) :: me logical :: ok ok = ( me%iflag == 0_ip ) end function status_ok !***************************************************************************************** !***************************************************************************************** !> ! This sets the `iflag` variable in the class to `0` ! (which indicates that everything is OK). It can be used ! after an error is encountered. elemental subroutine clear_bspline_flag(me) implicit none class(bspline_class),intent(inout) :: me me%iflag = 0_ip end subroutine clear_bspline_flag !***************************************************************************************** !***************************************************************************************** !> ! Get the status message from a [[bspline_class]] routine call. ! ! If `iflag` is not included, then the one in the class is used (which ! corresponds to the last routine called.) ! Otherwise, it will convert the ! input `iflag` argument into the appropriate message. ! ! This is a wrapper for [[get_status_message]]. pure function get_bspline_status_message(me,iflag) result(msg) implicit none class(bspline_class),intent(in) :: me character(len=:),allocatable :: msg !! status message associated with the flag integer(ip),intent(in),optional :: iflag !! the corresponding status code if (present(iflag)) then msg = get_status_message(iflag) else msg = get_status_message(me%iflag) end if end function get_bspline_status_message !***************************************************************************************** !***************************************************************************************** !> ! Actual size of a [[bspline_1d]] structure in bits. pure function size_1d(me) result(s) implicit none class(bspline_1d),intent(in) :: me integer(ip) :: s !! size of the structure in bits s = 3_ip*int_size + logical_size + 2_ip*int_size if (allocated(me%bcoef)) s = s + real_size*size(me%bcoef,kind=ip) if (allocated(me%tx)) s = s + real_size*size(me%tx,kind=ip) if (allocated(me%work_val_1)) s = s + real_size*size(me%work_val_1,kind=ip) end function size_1d !***************************************************************************************** !***************************************************************************************** !> ! Actual size of a [[bspline_2d]] structure in bits. pure function size_2d(me) result(s) implicit none class(bspline_2d),intent(in) :: me integer(ip) :: s !! size of the structure in bits s = 3_ip*int_size + logical_size + 6_ip*int_size if (allocated(me%bcoef)) s = s + real_size*size(me%bcoef,kind=ip) if (allocated(me%tx)) s = s + real_size*size(me%tx,kind=ip) if (allocated(me%ty)) s = s + real_size*size(me%ty,kind=ip) if (allocated(me%work_val_1)) s = s + real_size*size(me%work_val_1,kind=ip) if (allocated(me%work_val_2)) s = s + real_size*size(me%work_val_2,kind=ip) end function size_2d !***************************************************************************************** !***************************************************************************************** !> ! Actual size of a [[bspline_3d]] structure in bits. pure function size_3d(me) result(s) implicit none class(bspline_3d),intent(in) :: me integer(ip) :: s !! size of the structure in bits s = 3_ip*int_size + logical_size + 10_ip*int_size if (allocated(me%bcoef)) s = s + real_size*size(me%bcoef,kind=ip) if (allocated(me%tx)) s = s + real_size*size(me%tx,kind=ip) if (allocated(me%ty)) s = s + real_size*size(me%ty,kind=ip) if (allocated(me%tz)) s = s + real_size*size(me%tz,kind=ip) if (allocated(me%work_val_1)) s = s + real_size*size(me%work_val_1,1_ip,kind=ip)*& size(me%work_val_1,2_ip,kind=ip) if (allocated(me%work_val_2)) s = s + real_size*size(me%work_val_2,kind=ip) if (allocated(me%work_val_3)) s = s + real_size*size(me%work_val_3,kind=ip) end function size_3d !***************************************************************************************** !***************************************************************************************** !> ! Actual size of a [[bspline_4d]] structure in bits. pure function size_4d(me) result(s) implicit none class(bspline_4d),intent(in) :: me integer(ip) :: s !! size of the structure in bits s = 3_ip*int_size + logical_size + 14_ip*int_size if (allocated(me%bcoef)) s = s + real_size*size(me%bcoef,kind=ip) if (allocated(me%tx)) s = s + real_size*size(me%tx,kind=ip) if (allocated(me%ty)) s = s + real_size*size(me%ty,kind=ip) if (allocated(me%tz)) s = s + real_size*size(me%tz,kind=ip) if (allocated(me%tq)) s = s + real_size*size(me%tq,kind=ip) if (allocated(me%work_val_1)) s = s + real_size*size(me%work_val_1,1_ip,kind=ip)*& size(me%work_val_1,2_ip,kind=ip)*& size(me%work_val_1,3_ip,kind=ip) if (allocated(me%work_val_2)) s = s + real_size*size(me%work_val_2,1_ip,kind=ip)*& size(me%work_val_2,2_ip,kind=ip) if (allocated(me%work_val_3)) s = s + real_size*size(me%work_val_3,kind=ip) if (allocated(me%work_val_4)) s = s + real_size*size(me%work_val_4,kind=ip) end function size_4d !***************************************************************************************** !***************************************************************************************** !> ! Actual size of a [[bspline_5d]] structure in bits. pure function size_5d(me) result(s) implicit none class(bspline_5d),intent(in) :: me integer(ip) :: s !! size of the structure in bits s = 3_ip*int_size + logical_size + 18_ip*int_size if (allocated(me%bcoef)) s = s + real_size*size(me%bcoef,kind=ip) if (allocated(me%tx)) s = s + real_size*size(me%tx,kind=ip) if (allocated(me%ty)) s = s + real_size*size(me%ty,kind=ip) if (allocated(me%tz)) s = s + real_size*size(me%tz,kind=ip) if (allocated(me%tq)) s = s + real_size*size(me%tq,kind=ip) if (allocated(me%tr)) s = s + real_size*size(me%tr,kind=ip) if (allocated(me%work_val_1)) s = s + real_size*size(me%work_val_1,1_ip,kind=ip)*& size(me%work_val_1,2_ip,kind=ip)*& size(me%work_val_1,3_ip,kind=ip)*& size(me%work_val_1,4_ip,kind=ip) if (allocated(me%work_val_2)) s = s + real_size*size(me%work_val_2,1_ip,kind=ip)*& size(me%work_val_2,2_ip,kind=ip)*& size(me%work_val_2,3_ip,kind=ip) if (allocated(me%work_val_3)) s = s + real_size*size(me%work_val_3,1_ip,kind=ip)*& size(me%work_val_3,2_ip,kind=ip) if (allocated(me%work_val_4)) s = s + real_size*size(me%work_val_4,kind=ip) if (allocated(me%work_val_5)) s = s + real_size*size(me%work_val_5,kind=ip) end function size_5d !***************************************************************************************** !***************************************************************************************** !> ! Actual size of a [[bspline_6d]] structure in bits. pure function size_6d(me) result(s) implicit none class(bspline_6d),intent(in) :: me integer(ip) :: s !! size of the structure in bits s = 3_ip*int_size + logical_size + 22_ip*int_size if (allocated(me%bcoef)) s = s + real_size*size(me%bcoef,kind=ip) if (allocated(me%tx)) s = s + real_size*size(me%tx,kind=ip) if (allocated(me%ty)) s = s + real_size*size(me%ty,kind=ip) if (allocated(me%tz)) s = s + real_size*size(me%tz,kind=ip) if (allocated(me%tq)) s = s + real_size*size(me%tq,kind=ip) if (allocated(me%tr)) s = s + real_size*size(me%tr,kind=ip) if (allocated(me%ts)) s = s + real_size*size(me%ts,kind=ip) if (allocated(me%work_val_1)) s = s + real_size*size(me%work_val_1,1_ip,kind=ip)*& size(me%work_val_1,2_ip,kind=ip)*& size(me%work_val_1,3_ip,kind=ip)*& size(me%work_val_1,4_ip,kind=ip)*& size(me%work_val_1,5_ip,kind=ip) if (allocated(me%work_val_2)) s = s + real_size*size(me%work_val_2,1_ip,kind=ip)*& size(me%work_val_2,2_ip,kind=ip)*& size(me%work_val_2,3_ip,kind=ip)*& size(me%work_val_2,4_ip,kind=ip) if (allocated(me%work_val_3)) s = s + real_size*size(me%work_val_3,1_ip,kind=ip)*& size(me%work_val_3,2_ip,kind=ip)*& size(me%work_val_3,3_ip,kind=ip) if (allocated(me%work_val_4)) s = s + real_size*size(me%work_val_4,1_ip,kind=ip)*& size(me%work_val_4,2_ip,kind=ip) if (allocated(me%work_val_5)) s = s + real_size*size(me%work_val_5,kind=ip) if (allocated(me%work_val_6)) s = s + real_size*size(me%work_val_6,kind=ip) end function size_6d !***************************************************************************************** !***************************************************************************************** !> ! Destructor for contents of the base [[bspline_class]] class. ! (this routine is called by the extended classes). pure subroutine destroy_base(me) implicit none class(bspline_class),intent(inout) :: me me%inbvx = 1_ip me%iflag = 1_ip me%initialized = .false. me%extrap = .false. me%nf = 0_ip end subroutine destroy_base !***************************************************************************************** !***************************************************************************************** !> ! Destructor for [[bspline_1d]] class. pure subroutine destroy_1d(me) implicit none class(bspline_1d),intent(inout) :: me call me%destroy_base() me%nx = 0_ip me%kx = 0_ip if (allocated(me%bcoef)) deallocate(me%bcoef) if (allocated(me%tx)) deallocate(me%tx) if (allocated(me%work_val_1)) deallocate(me%work_val_1) end subroutine destroy_1d !***************************************************************************************** !***************************************************************************************** !> ! Destructor for [[bspline_2d]] class. pure subroutine destroy_2d(me) implicit none class(bspline_2d),intent(inout) :: me call me%destroy_base() me%nx = 0_ip me%ny = 0_ip me%kx = 0_ip me%ky = 0_ip me%inbvy = 1_ip me%iloy = 1_ip if (allocated(me%bcoef)) deallocate(me%bcoef) if (allocated(me%tx)) deallocate(me%tx) if (allocated(me%ty)) deallocate(me%ty) if (allocated(me%work_val_1)) deallocate(me%work_val_1) if (allocated(me%work_val_2)) deallocate(me%work_val_2) end subroutine destroy_2d !***************************************************************************************** !***************************************************************************************** !> ! Destructor for [[bspline_3d]] class. pure subroutine destroy_3d(me) implicit none class(bspline_3d),intent(inout) :: me call me%destroy_base() me%nx = 0_ip me%ny = 0_ip me%nz = 0_ip me%kx = 0_ip me%ky = 0_ip me%kz = 0_ip me%inbvy = 1_ip me%inbvz = 1_ip me%iloy = 1_ip me%iloz = 1_ip if (allocated(me%bcoef)) deallocate(me%bcoef) if (allocated(me%tx)) deallocate(me%tx) if (allocated(me%ty)) deallocate(me%ty) if (allocated(me%tz)) deallocate(me%tz) if (allocated(me%work_val_1)) deallocate(me%work_val_1) if (allocated(me%work_val_2)) deallocate(me%work_val_2) if (allocated(me%work_val_3)) deallocate(me%work_val_3) end subroutine destroy_3d !***************************************************************************************** !***************************************************************************************** !> ! Destructor for [[bspline_4d]] class. pure subroutine destroy_4d(me) implicit none class(bspline_4d),intent(inout) :: me call me%destroy_base() me%nx = 0_ip me%ny = 0_ip me%nz = 0_ip me%nq = 0_ip me%kx = 0_ip me%ky = 0_ip me%kz = 0_ip me%kq = 0_ip me%inbvy = 1_ip me%inbvz = 1_ip me%inbvq = 1_ip me%iloy = 1_ip me%iloz = 1_ip me%iloq = 1_ip if (allocated(me%bcoef)) deallocate(me%bcoef) if (allocated(me%tx)) deallocate(me%tx) if (allocated(me%ty)) deallocate(me%ty) if (allocated(me%tz)) deallocate(me%tz) if (allocated(me%tq)) deallocate(me%tq) if (allocated(me%work_val_1)) deallocate(me%work_val_1) if (allocated(me%work_val_2)) deallocate(me%work_val_2) if (allocated(me%work_val_3)) deallocate(me%work_val_3) if (allocated(me%work_val_4)) deallocate(me%work_val_4) end subroutine destroy_4d !***************************************************************************************** !***************************************************************************************** !> ! Destructor for [[bspline_5d]] class. pure subroutine destroy_5d(me) implicit none class(bspline_5d),intent(inout) :: me call me%destroy_base() me%nx = 0_ip me%ny = 0_ip me%nz = 0_ip me%nq = 0_ip me%nr = 0_ip me%kx = 0_ip me%ky = 0_ip me%kz = 0_ip me%kq = 0_ip me%kr = 0_ip me%inbvy = 1_ip me%inbvz = 1_ip me%inbvq = 1_ip me%inbvr = 1_ip me%iloy = 1_ip me%iloz = 1_ip me%iloq = 1_ip me%ilor = 1_ip if (allocated(me%bcoef)) deallocate(me%bcoef) if (allocated(me%tx)) deallocate(me%tx) if (allocated(me%ty)) deallocate(me%ty) if (allocated(me%tz)) deallocate(me%tz) if (allocated(me%tq)) deallocate(me%tq) if (allocated(me%tr)) deallocate(me%tr) if (allocated(me%work_val_1)) deallocate(me%work_val_1) if (allocated(me%work_val_2)) deallocate(me%work_val_2) if (allocated(me%work_val_3)) deallocate(me%work_val_3) if (allocated(me%work_val_4)) deallocate(me%work_val_4) if (allocated(me%work_val_5)) deallocate(me%work_val_5) end subroutine destroy_5d !***************************************************************************************** !***************************************************************************************** !> ! Destructor for [[bspline_6d]] class. pure subroutine destroy_6d(me) implicit none class(bspline_6d),intent(inout) :: me call me%destroy_base() me%nx = 0_ip me%ny = 0_ip me%nz = 0_ip me%nq = 0_ip me%nr = 0_ip me%ns = 0_ip me%kx = 0_ip me%ky = 0_ip me%kz = 0_ip me%kq = 0_ip me%kr = 0_ip me%ks = 0_ip me%inbvy = 1_ip me%inbvz = 1_ip me%inbvq = 1_ip me%inbvr = 1_ip me%inbvs = 1_ip me%iloy = 1_ip me%iloz = 1_ip me%iloq = 1_ip me%ilor = 1_ip me%ilos = 1_ip if (allocated(me%bcoef)) deallocate(me%bcoef) if (allocated(me%tx)) deallocate(me%tx) if (allocated(me%ty)) deallocate(me%ty) if (allocated(me%tz)) deallocate(me%tz) if (allocated(me%tq)) deallocate(me%tq) if (allocated(me%tr)) deallocate(me%tr) if (allocated(me%ts)) deallocate(me%ts) if (allocated(me%work_val_1)) deallocate(me%work_val_1) if (allocated(me%work_val_2)) deallocate(me%work_val_2) if (allocated(me%work_val_3)) deallocate(me%work_val_3) if (allocated(me%work_val_4)) deallocate(me%work_val_4) if (allocated(me%work_val_5)) deallocate(me%work_val_5) if (allocated(me%work_val_6)) deallocate(me%work_val_6) end subroutine destroy_6d !***************************************************************************************** !***************************************************************************************** !> ! Finalizer for [[bspline_1d]] class. Just a wrapper for [[destroy_1d]]. pure elemental subroutine finalize_1d(me) type(bspline_1d),intent(inout) :: me; call me%destroy() end subroutine finalize_1d !***************************************************************************************** !***************************************************************************************** !> ! Finalizer for [[bspline_2d]] class. Just a wrapper for [[destroy_2d]]. pure elemental subroutine finalize_2d(me) type(bspline_2d),intent(inout) :: me; call me%destroy() end subroutine finalize_2d !***************************************************************************************** !***************************************************************************************** !> ! Finalizer for [[bspline_3d]] class. Just a wrapper for [[destroy_3d]]. pure elemental subroutine finalize_3d(me) type(bspline_3d),intent(inout) :: me; call me%destroy() end subroutine finalize_3d !***************************************************************************************** !***************************************************************************************** !> ! Finalizer for [[bspline_4d]] class. Just a wrapper for [[destroy_4d]]. pure elemental subroutine finalize_4d(me) type(bspline_4d),intent(inout) :: me; call me%destroy() end subroutine finalize_4d !***************************************************************************************** !***************************************************************************************** !> ! Finalizer for [[bspline_5d]] class. Just a wrapper for [[destroy_5d]]. pure elemental subroutine finalize_5d(me) type(bspline_5d),intent(inout) :: me; call me%destroy() end subroutine finalize_5d !***************************************************************************************** !***************************************************************************************** !> ! Finalizer for [[bspline_6d]] class. Just a wrapper for [[destroy_6d]]. pure elemental subroutine finalize_6d(me) type(bspline_6d),intent(inout) :: me; call me%destroy() end subroutine finalize_6d !***************************************************************************************** !***************************************************************************************** !> ! Sets the `extrap` flag in the class. pure subroutine set_extrap_flag(me,extrap) implicit none class(bspline_class),intent(inout) :: me logical,intent(in),optional :: extrap !! if not present, then False is used if (present(extrap)) then me%extrap = extrap else me%extrap = .false. end if end subroutine set_extrap_flag !***************************************************************************************** !***************************************************************************************** !> ! It returns an empty [[bspline_1d]] type. Note that INITIALIZE still ! needs to be called before it can be used. ! Not really that useful except perhaps in some OpenMP applications. pure elemental function bspline_1d_constructor_empty() result(me) implicit none type(bspline_1d) :: me end function bspline_1d_constructor_empty !***************************************************************************************** !***************************************************************************************** !> ! Constructor for a [[bspline_1d]] type (auto knots). ! This is a wrapper for [[initialize_1d_auto_knots]]. pure function bspline_1d_constructor_auto_knots(x,fcn,kx,extrap) result(me) implicit none type(bspline_1d) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: fcn !! `(nx)` array of function values to interpolate. `fcn(i)` should !! contain the function value at the point `x(i)` integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) call initialize_1d_auto_knots(me,x,fcn,kx,me%iflag,extrap) end function bspline_1d_constructor_auto_knots !***************************************************************************************** !***************************************************************************************** !> ! Constructor for a [[bspline_1d]] type (user-specified knots). ! This is a wrapper for [[initialize_1d_specify_knots]]. pure function bspline_1d_constructor_specify_knots(x,fcn,kx,tx,extrap) result(me) implicit none type(bspline_1d) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: fcn !! `(nx)` array of function values to interpolate. `fcn(i)` should !! contain the function value at the point `x(i)` integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction !! for the spline interpolant. !! Must be non-decreasing. logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) call initialize_1d_specify_knots(me,x,fcn,kx,tx,me%iflag,extrap) end function bspline_1d_constructor_specify_knots !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_1d]] type (with automatically-computed knots). ! This is a wrapper for [[db1ink]]. pure subroutine initialize_1d_auto_knots(me,x,fcn,kx,iflag,extrap) implicit none class(bspline_1d),intent(inout) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: fcn !! `(nx)` array of function values to interpolate. `fcn(i)` should !! contain the function value at the point `x(i)` integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(out) :: iflag !! status flag (see [[db1ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip) :: iknot integer(ip) :: nx call me%destroy() nx = size(x,kind=ip) me%nx = nx me%kx = kx allocate(me%tx(nx+kx)) allocate(me%bcoef(nx,1_ip)) allocate(me%work_val_1(3_ip*kx)) iknot = 0_ip !knot sequence chosen by db1ink call db1ink(x,nx,fcn,kx,iknot,me%tx,me%bcoef(:,1),iflag) if (iflag==0_ip) then call me%set_extrap_flag(extrap) end if if (iflag==0_ip) me%nf = 1_ip me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_1d_auto_knots !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_1d]] type (with user-specified knots). ! This is a wrapper for [[db1ink]]. pure subroutine initialize_1d_specify_knots(me,x,fcn,kx,tx,iflag,extrap) implicit none class(bspline_1d),intent(inout) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: fcn !! `(nx)` array of function values to interpolate. `fcn(i)` should !! contain the function value at the point `x(i)` integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction !! for the spline interpolant. !! Must be non-decreasing. integer(ip),intent(out) :: iflag !! status flag (see [[db1ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip) :: nx call me%destroy() nx = size(x,kind=ip) call check_knot_vectors_sizes(nx=nx,kx=kx,tx=tx,iflag=iflag) if (iflag == 0_ip) then me%nx = nx me%kx = kx allocate(me%tx(nx+kx)) allocate(me%bcoef(nx,1_ip)) allocate(me%work_val_1(3_ip*kx)) me%tx = tx call db1ink(x,nx,fcn,kx,1_ip,me%tx,me%bcoef(:,1),iflag) call me%set_extrap_flag(extrap) end if if (iflag==0_ip) me%nf = 1_ip me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_1d_specify_knots !***************************************************************************************** !***************************************************************************************** !> ! Evaluate a [[bspline_1d]] interpolate. This is a wrapper for [[db1val]]. pure subroutine evaluate_1d(me,xval,idx,f,iflag) implicit none class(bspline_1d),intent(inout) :: me real(wp),intent(in) :: xval !! \(x\) coordinate of evaluation point. integer(ip),intent(in) :: idx !! \(x\) derivative of piecewise polynomial to evaluate. real(wp),intent(out) :: f !! interpolated value integer(ip),intent(out) :: iflag !! status flag (see [[db1val]]) if (me%initialized .and. me%nf/=1_ip) then f = 0.0_wp iflag = 6003_ip ! use the multiple-output version else if (me%initialized) then call db1val(xval,idx,me%tx,me%nx,me%kx,me%bcoef(:,1),f,iflag,& me%inbvx,me%work_val_1,extrap=me%extrap) else f = 0.0_wp iflag = 1_ip end if me%iflag = iflag end subroutine evaluate_1d !***************************************************************************************** !***************************************************************************************** !> ! Evaluate a [[bspline_1d]] definite integral. This is a wrapper for [[db1sqad]]. pure subroutine integral_1d(me,x1,x2,f,iflag) implicit none class(bspline_1d),intent(inout) :: me real(wp),intent(in) :: x1 !! left point of interval real(wp),intent(in) :: x2 !! right point of interval real(wp),intent(out) :: f !! integral of the b-spline over \( [x_1, x_2] \) integer(ip),intent(out) :: iflag !! status flag (see [[db1sqad]]) if (me%initialized .and. me%nf/=1_ip) then iflag = 6003_ip ! use the multiple-output version else if (me%initialized) then call db1sqad(me%tx,me%bcoef(:,1),me%nx,me%kx,x1,x2,f,iflag,me%work_val_1) else iflag = 1_ip end if me%iflag = iflag end subroutine integral_1d !***************************************************************************************** !***************************************************************************************** !> ! Evaluate a [[bspline_1d]] definite integral. This is a wrapper for [[db1fqad]]. subroutine fintegral_1d(me,fun,idx,x1,x2,tol,f,iflag) implicit none class(bspline_1d),intent(inout) :: me procedure(b1fqad_func) :: fun !! external function of one argument for the !! integrand `bf(x)=fun(x)*dbvalu(tx,bcoef,nx,kx,idx,x,inbv)` integer(ip),intent(in) :: idx !! order of the spline derivative, `0 <= idx <= k-1` !! `idx=0` gives the spline function real(wp),intent(in) :: x1 !! left point of interval real(wp),intent(in) :: x2 !! right point of interval real(wp),intent(in) :: tol !! desired accuracy for the quadrature real(wp),intent(out) :: f !! integral of `bf(x)` over \( [x_1, x_2] \) integer(ip),intent(out) :: iflag !! status flag (see [[db1sqad]]) if (me%initialized .and. me%nf/=1_ip) then iflag = 6003_ip ! use the multiple-output version else if (me%initialized) then call db1fqad(fun,me%tx,me%bcoef(:,1),me%nx,me%kx,idx,x1,x2,tol,f,iflag,me%work_val_1) else iflag = 1_ip end if me%iflag = iflag end subroutine fintegral_1d !***************************************************************************************** !***************************************************************************************** !> ! All the roots of \( s(x) = a \) for a [[bspline_1d]]. ! This is a wrapper for [[db1root]]. pure subroutine roots_1d(me,a,x,iflag) implicit none class(bspline_1d),intent(inout) :: me !! the spline real(wp),intent(in) :: a !! the value to solve for real(wp),dimension(:),allocatable,intent(out) :: x !! the roots, in increasing order !! (size 0 if there are none) integer(ip),intent(out) :: iflag !! status flag (see [[db1root]]) if (me%initialized .and. me%nf/=1_ip) then iflag = 6003_ip ! use the multiple-output version allocate(x(0)) else if (me%initialized) then call db1root(me%tx,me%nx,me%kx,me%bcoef(:,1),a,x,iflag) else allocate(x(0)) iflag = 1_ip end if me%iflag = iflag end subroutine roots_1d !***************************************************************************************** !***************************************************************************************** !> ! The unique \( x \) with \( s(x) = y \) for a [[bspline_1d]] ! (for example, a monotonic spline). This is a wrapper for [[db1root]]. ! ! Returns `iflag=4005` if there is no solution, and `iflag=4006` ! if there is more than one (in which case `x` is the smallest). pure subroutine inverse_1d(me,y,x,iflag) implicit none class(bspline_1d),intent(inout) :: me !! the spline real(wp),intent(in) :: y !! the spline value real(wp),intent(out) :: x !! the solution of \( s(x) = y \) integer(ip),intent(out) :: iflag !! status flag (see [[db1root]]) real(wp),dimension(:),allocatable :: roots !! all the roots of \( s(x) = y \) x = 0.0_wp call me%roots(y,roots,iflag) if (iflag==0_ip .or. iflag==4004_ip) then if (size(roots)>0) x = roots(1) if (iflag==4004_ip .or. size(roots)>1) then iflag = 4006_ip ! s(x) = y on a whole interval, or more than one root else if (size(roots)==0) then iflag = 4005_ip end if end if me%iflag = iflag end subroutine inverse_1d !***************************************************************************************** !***************************************************************************************** !> ! The global minimum of a [[bspline_1d]] on \( [x_a, x_b] \) ! (the whole spline domain if `xa` and `xb` are not present). ! This is a wrapper for [[db1min]]. pure subroutine minimum_1d(me,xmin,fmin,iflag,xa,xb) implicit none class(bspline_1d),intent(inout) :: me !! the spline real(wp),intent(out) :: xmin !! location of the minimum real(wp),intent(out) :: fmin !! minimum value integer(ip),intent(out) :: iflag !! status flag (see [[db1min]]) real(wp),intent(in),optional :: xa !! left end of the interval (default `tx(kx)`) real(wp),intent(in),optional :: xb !! right end of the interval (default `tx(nx+1)`) real(wp) :: x1 !! left end of the interval real(wp) :: x2 !! right end of the interval if (me%initialized .and. me%nf/=1_ip) then iflag = 6003_ip ! use the multiple-output version xmin = 0.0_wp fmin = 0.0_wp else if (me%initialized) then call get_interval_1d(me,xa,xb,x1,x2) call db1min(me%tx,me%nx,me%kx,me%bcoef(:,1),x1,x2,xmin,fmin,iflag) else xmin = 0.0_wp fmin = 0.0_wp iflag = 1_ip end if me%iflag = iflag end subroutine minimum_1d !***************************************************************************************** !***************************************************************************************** !> ! The global maximum of a [[bspline_1d]] on \( [x_a, x_b] \) ! (the whole spline domain if `xa` and `xb` are not present). ! This is a wrapper for [[db1max]]. pure subroutine maximum_1d(me,xmax,fmax,iflag,xa,xb) implicit none class(bspline_1d),intent(inout) :: me !! the spline real(wp),intent(out) :: xmax !! location of the maximum real(wp),intent(out) :: fmax !! maximum value integer(ip),intent(out) :: iflag !! status flag (see [[db1max]]) real(wp),intent(in),optional :: xa !! left end of the interval (default `tx(kx)`) real(wp),intent(in),optional :: xb !! right end of the interval (default `tx(nx+1)`) real(wp) :: x1 !! left end of the interval real(wp) :: x2 !! right end of the interval if (me%initialized .and. me%nf/=1_ip) then iflag = 6003_ip ! use the multiple-output version xmax = 0.0_wp fmax = 0.0_wp else if (me%initialized) then call get_interval_1d(me,xa,xb,x1,x2) call db1max(me%tx,me%nx,me%kx,me%bcoef(:,1),x1,x2,xmax,fmax,iflag) else xmax = 0.0_wp fmax = 0.0_wp iflag = 1_ip end if me%iflag = iflag end subroutine maximum_1d !***************************************************************************************** !***************************************************************************************** !> ! The interval for [[minimum_1d]] and [[maximum_1d]] ! (the whole spline domain by default). pure subroutine get_interval_1d(me,xa,xb,x1,x2) implicit none class(bspline_1d),intent(in) :: me !! the spline real(wp),intent(in),optional :: xa !! left end of the interval real(wp),intent(in),optional :: xb !! right end of the interval real(wp),intent(out) :: x1 !! `xa`, or `tx(kx)` if not present real(wp),intent(out) :: x2 !! `xb`, or `tx(nx+1)` if not present if (present(xa)) then x1 = xa else x1 = me%tx(me%kx) end if if (present(xb)) then x2 = xb else x2 = me%tx(me%nx+1_ip) end if end subroutine get_interval_1d !***************************************************************************************** !***************************************************************************************** !> ! It returns an empty [[bspline_2d]] type. Note that INITIALIZE still ! needs to be called before it can be used. ! Not really that useful except perhaps in some OpenMP applications. elemental function bspline_2d_constructor_empty() result(me) implicit none type(bspline_2d) :: me end function bspline_2d_constructor_empty !***************************************************************************************** !***************************************************************************************** !> ! Constructor for a [[bspline_2d]] type (auto knots). ! This is a wrapper for [[initialize_2d_auto_knots]]. pure function bspline_2d_constructor_auto_knots(x,y,fcn,kx,ky,extrap) result(me) implicit none type(bspline_2d) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:,:),intent(in) :: fcn !! `(nx,ny)` matrix of function values to interpolate. !! `fcn(i,j)` should contain the function value at the !! point (`x(i)`,`y(j)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) call initialize_2d_auto_knots(me,x,y,fcn,kx,ky,me%iflag,extrap) end function bspline_2d_constructor_auto_knots !***************************************************************************************** !***************************************************************************************** !> ! Constructor for a [[bspline_2d]] type (user-specified knots). ! This is a wrapper for [[initialize_2d_specify_knots]]. pure function bspline_2d_constructor_specify_knots(x,y,fcn,kx,ky,tx,ty,extrap) result(me) implicit none type(bspline_2d) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:,:),intent(in) :: fcn !! `(nx,ny)` matrix of function values to interpolate. !! `fcn(i,j)` should contain the function value at the !! point (`x(i)`,`y(j)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction !! for the spline interpolant. !! Must be non-decreasing. logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) call initialize_2d_specify_knots(me,x,y,fcn,kx,ky,tx,ty,me%iflag,extrap) end function bspline_2d_constructor_specify_knots !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_2d]] type (with automatically-computed knots). ! This is a wrapper for [[db2ink]]. pure subroutine initialize_2d_auto_knots(me,x,y,fcn,kx,ky,iflag,extrap) implicit none class(bspline_2d),intent(inout) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:,:),intent(in) :: fcn !! `(nx,ny)` matrix of function values to interpolate. !! `fcn(i,j)` should contain the function value at the !! point (`x(i)`,`y(j)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(out) :: iflag !! status flag (see [[db2ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip) :: iknot integer(ip) :: nx,ny call me%destroy() nx = size(x,kind=ip) ny = size(y,kind=ip) me%nx = nx me%ny = ny me%kx = kx me%ky = ky allocate(me%tx(nx+kx)) allocate(me%ty(ny+ky)) allocate(me%bcoef(nx,ny,1_ip)) allocate(me%work_val_1(ky)) allocate(me%work_val_2(3_ip*max(kx,ky))) iknot = 0_ip !knot sequence chosen by db2ink call db2ink(x,nx,y,ny,fcn,kx,ky,iknot,me%tx,me%ty,me%bcoef(:,:,1),iflag) if (iflag==0_ip) then call me%set_extrap_flag(extrap) end if if (iflag==0_ip) me%nf = 1_ip me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_2d_auto_knots !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_2d]] type (with user-specified knots). ! This is a wrapper for [[db2ink]]. pure subroutine initialize_2d_specify_knots(me,x,y,fcn,kx,ky,tx,ty,iflag,extrap) implicit none class(bspline_2d),intent(inout) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:,:),intent(in) :: fcn !! `(nx,ny)` matrix of function values to interpolate. !! `fcn(i,j)` should contain the function value at the !! point (`x(i)`,`y(j)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction !! for the spline interpolant. !! Must be non-decreasing. integer(ip),intent(out) :: iflag !! status flag (see [[db2ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip) :: nx,ny call me%destroy() nx = size(x,kind=ip) ny = size(y,kind=ip) call check_knot_vectors_sizes(nx=nx,kx=kx,tx=tx,& ny=ny,ky=ky,ty=ty,& iflag=iflag) if (iflag == 0_ip) then me%nx = nx me%ny = ny me%kx = kx me%ky = ky allocate(me%tx(nx+kx)) allocate(me%ty(ny+ky)) allocate(me%bcoef(nx,ny,1_ip)) allocate(me%work_val_1(ky)) allocate(me%work_val_2(3_ip*max(kx,ky))) me%tx = tx me%ty = ty call db2ink(x,nx,y,ny,fcn,kx,ky,1_ip,me%tx,me%ty,me%bcoef(:,:,1),iflag) call me%set_extrap_flag(extrap) end if if (iflag==0_ip) me%nf = 1_ip me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_2d_specify_knots !***************************************************************************************** !***************************************************************************************** !> ! Evaluate a [[bspline_2d]] interpolate. This is a wrapper for [[db2val]]. pure subroutine evaluate_2d(me,xval,yval,idx,idy,f,iflag) implicit none class(bspline_2d),intent(inout) :: me real(wp),intent(in) :: xval !! \(x\) coordinate of evaluation point. real(wp),intent(in) :: yval !! \(y\) coordinate of evaluation point. integer(ip),intent(in) :: idx !! \(x\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idy !! \(y\) derivative of piecewise polynomial to evaluate. real(wp),intent(out) :: f !! interpolated value integer(ip),intent(out) :: iflag !! status flag (see [[db2val]]) if (me%initialized .and. me%nf/=1_ip) then f = 0.0_wp iflag = 6003_ip ! use the multiple-output version else if (me%initialized) then call db2val(xval,yval,& idx,idy,& me%tx,me%ty,& me%nx,me%ny,& me%kx,me%ky,& me%bcoef(:,:,1),f,iflag,& me%inbvx,me%inbvy,me%iloy,& me%work_val_1,me%work_val_2,& extrap=me%extrap) else f = 0.0_wp iflag = 1_ip end if me%iflag = iflag end subroutine evaluate_2d !***************************************************************************************** !***************************************************************************************** !> ! It returns an empty [[bspline_3d]] type. Note that INITIALIZE still ! needs to be called before it can be used. ! Not really that useful except perhaps in some OpenMP applications. elemental function bspline_3d_constructor_empty() result(me) implicit none type(bspline_3d) :: me end function bspline_3d_constructor_empty !***************************************************************************************** !***************************************************************************************** !> ! Constructor for a [[bspline_3d]] type (auto knots). ! This is a wrapper for [[initialize_3d_auto_knots]]. pure function bspline_3d_constructor_auto_knots(x,y,z,fcn,kx,ky,kz,extrap) result(me) implicit none type(bspline_3d) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:),intent(in) :: fcn !! `(nx,ny,nz)` matrix of function values to interpolate. !! `fcn(i,j,k)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) call initialize_3d_auto_knots(me,x,y,z,fcn,kx,ky,kz,me%iflag,extrap) end function bspline_3d_constructor_auto_knots !***************************************************************************************** !***************************************************************************************** !> ! Constructor for a [[bspline_3d]] type (user-specified knots). ! This is a wrapper for [[initialize_3d_specify_knots]]. pure function bspline_3d_constructor_specify_knots(x,y,z,fcn,kx,ky,kz,tx,ty,tz,extrap) result(me) implicit none type(bspline_3d) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:),intent(in) :: fcn !! `(nx,ny,nz)` matrix of function values to interpolate. !! `fcn(i,j,k)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction !! for the spline interpolant. !! Must be non-decreasing. logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) call initialize_3d_specify_knots(me,x,y,z,fcn,kx,ky,kz,tx,ty,tz,me%iflag,extrap) end function bspline_3d_constructor_specify_knots !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_3d]] type (with automatically-computed knots). ! This is a wrapper for [[db3ink]]. pure subroutine initialize_3d_auto_knots(me,x,y,z,fcn,kx,ky,kz,iflag,extrap) implicit none class(bspline_3d),intent(inout) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:),intent(in) :: fcn !! `(nx,ny,nz)` matrix of function values to interpolate. !! `fcn(i,j,k)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(out) :: iflag !! status flag (see [[db3ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip) :: iknot integer(ip) :: nx,ny,nz call me%destroy() nx = size(x,kind=ip) ny = size(y,kind=ip) nz = size(z,kind=ip) me%nx = nx me%ny = ny me%nz = nz me%kx = kx me%ky = ky me%kz = kz allocate(me%tx(nx+kx)) allocate(me%ty(ny+ky)) allocate(me%tz(nz+kz)) allocate(me%bcoef(nx,ny,nz,1_ip)) allocate(me%work_val_1(ky,kz)) allocate(me%work_val_2(kz)) allocate(me%work_val_3(3_ip*max(kx,ky,kz))) iknot = 0_ip !knot sequence chosen by db3ink call db3ink(x,nx,y,ny,z,nz,& fcn,& kx,ky,kz,& iknot,& me%tx,me%ty,me%tz,& me%bcoef(:,:,:,1),iflag) if (iflag==0_ip) then call me%set_extrap_flag(extrap) end if if (iflag==0_ip) me%nf = 1_ip me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_3d_auto_knots !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_3d]] type (with user-specified knots). ! This is a wrapper for [[db3ink]]. pure subroutine initialize_3d_specify_knots(me,x,y,z,fcn,kx,ky,kz,tx,ty,tz,iflag,extrap) implicit none class(bspline_3d),intent(inout) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:),intent(in) :: fcn !! `(nx,ny,nz)` matrix of function values to interpolate. !! `fcn(i,j,k)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction !! for the spline interpolant. !! Must be non-decreasing. integer(ip),intent(out) :: iflag !! status flag (see [[db3ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip) :: nx,ny,nz call me%destroy() nx = size(x,kind=ip) ny = size(y,kind=ip) nz = size(z,kind=ip) call check_knot_vectors_sizes(nx=nx,kx=kx,tx=tx,& ny=ny,ky=ky,ty=ty,& nz=nz,kz=kz,tz=tz,& iflag=iflag) if (iflag == 0_ip) then me%nx = nx me%ny = ny me%nz = nz me%kx = kx me%ky = ky me%kz = kz allocate(me%tx(nx+kx)) allocate(me%ty(ny+ky)) allocate(me%tz(nz+kz)) allocate(me%bcoef(nx,ny,nz,1_ip)) allocate(me%work_val_1(ky,kz)) allocate(me%work_val_2(kz)) allocate(me%work_val_3(3_ip*max(kx,ky,kz))) me%tx = tx me%ty = ty me%tz = tz call db3ink(x,nx,y,ny,z,nz,& fcn,& kx,ky,kz,& 1_ip,& me%tx,me%ty,me%tz,& me%bcoef(:,:,:,1),iflag) call me%set_extrap_flag(extrap) end if if (iflag==0_ip) me%nf = 1_ip me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_3d_specify_knots !***************************************************************************************** !***************************************************************************************** !> ! Evaluate a [[bspline_3d]] interpolate. This is a wrapper for [[db3val]]. pure subroutine evaluate_3d(me,xval,yval,zval,idx,idy,idz,f,iflag) implicit none class(bspline_3d),intent(inout) :: me real(wp),intent(in) :: xval !! \(x\) coordinate of evaluation point. real(wp),intent(in) :: yval !! \(y\) coordinate of evaluation point. real(wp),intent(in) :: zval !! \(z\) coordinate of evaluation point. integer(ip),intent(in) :: idx !! \(x\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idy !! \(y\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idz !! \(z\) derivative of piecewise polynomial to evaluate. real(wp),intent(out) :: f !! interpolated value integer(ip),intent(out) :: iflag !! status flag (see [[db3val]]) if (me%initialized .and. me%nf/=1_ip) then f = 0.0_wp iflag = 6003_ip ! use the multiple-output version else if (me%initialized) then call db3val(xval,yval,zval,& idx,idy,idz,& me%tx,me%ty,me%tz,& me%nx,me%ny,me%nz,& me%kx,me%ky,me%kz,& me%bcoef(:,:,:,1),f,iflag,& me%inbvx,me%inbvy,me%inbvz,& me%iloy,me%iloz,& me%work_val_1,me%work_val_2,me%work_val_3,& extrap=me%extrap) else f = 0.0_wp iflag = 1_ip end if me%iflag = iflag end subroutine evaluate_3d !***************************************************************************************** !***************************************************************************************** !> ! It returns an empty [[bspline_4d]] type. Note that INITIALIZE still ! needs to be called before it can be used. ! Not really that useful except perhaps in some OpenMP applications. elemental function bspline_4d_constructor_empty() result(me) implicit none type(bspline_4d) :: me end function bspline_4d_constructor_empty !***************************************************************************************** !***************************************************************************************** !> ! Constructor for a [[bspline_4d]] type (auto knots). ! This is a wrapper for [[initialize_4d_auto_knots]]. pure function bspline_4d_constructor_auto_knots(x,y,z,q,fcn,kx,ky,kz,kq,extrap) result(me) implicit none type(bspline_4d) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq)` matrix of function values to interpolate. !! `fcn(i,j,k,l)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 2 \le k_q < n_q \) ) !! (order = polynomial degree + 1) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) call initialize_4d_auto_knots(me,x,y,z,q,fcn,kx,ky,kz,kq,me%iflag,extrap) end function bspline_4d_constructor_auto_knots !***************************************************************************************** !***************************************************************************************** !> ! Constructor for a [[bspline_4d]] type (user-specified knots). ! This is a wrapper for [[initialize_4d_specify_knots]]. pure function bspline_4d_constructor_specify_knots(x,y,z,q,fcn,kx,ky,kz,kq,& tx,ty,tz,tq,extrap) result(me) implicit none type(bspline_4d) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq)` matrix of function values to interpolate. !! `fcn(i,j,k,l)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 2 \le k_q < n_q \) ) !! (order = polynomial degree + 1) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction !! for the spline interpolant. !! Must be non-decreasing. logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) call initialize_4d_specify_knots(me,x,y,z,q,fcn,kx,ky,kz,kq,tx,ty,tz,tq,me%iflag,extrap) end function bspline_4d_constructor_specify_knots !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_4d]] type (with automatically-computed knots). ! This is a wrapper for [[db4ink]]. pure subroutine initialize_4d_auto_knots(me,x,y,z,q,fcn,kx,ky,kz,kq,iflag,extrap) implicit none class(bspline_4d),intent(inout) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq)` matrix of function values to interpolate. !! `fcn(i,j,k,l)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 2 \le k_q < n_q \) ) !! (order = polynomial degree + 1) integer(ip),intent(out) :: iflag !! status flag (see [[db4ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip) :: iknot integer(ip) :: nx,ny,nz,nq call me%destroy() nx = size(x,kind=ip) ny = size(y,kind=ip) nz = size(z,kind=ip) nq = size(q,kind=ip) me%nx = nx me%ny = ny me%nz = nz me%nq = nq me%kx = kx me%ky = ky me%kz = kz me%kq = kq allocate(me%tx(nx+kx)) allocate(me%ty(ny+ky)) allocate(me%tz(nz+kz)) allocate(me%tq(nq+kq)) allocate(me%bcoef(nx,ny,nz,nq,1_ip)) allocate(me%work_val_1(ky,kz,kq)) allocate(me%work_val_2(kz,kq)) allocate(me%work_val_3(kq)) allocate(me%work_val_4(3_ip*max(kx,ky,kz,kq))) iknot = 0_ip !knot sequence chosen by db4ink call db4ink(x,nx,y,ny,z,nz,q,nq,& fcn,& kx,ky,kz,kq,& iknot,& me%tx,me%ty,me%tz,me%tq,& me%bcoef(:,:,:,:,1),iflag) if (iflag==0_ip) then call me%set_extrap_flag(extrap) end if if (iflag==0_ip) me%nf = 1_ip me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_4d_auto_knots !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_4d]] type (with user-specified knots). ! This is a wrapper for [[db4ink]]. pure subroutine initialize_4d_specify_knots(me,x,y,z,q,fcn,& kx,ky,kz,kq,tx,ty,tz,tq,iflag,extrap) implicit none class(bspline_4d),intent(inout) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq)` matrix of function values to interpolate. !! `fcn(i,j,k,l)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 2 \le k_q < n_q \) ) !! (order = polynomial degree + 1) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction !! for the spline interpolant. !! Must be non-decreasing. integer(ip),intent(out) :: iflag !! status flag (see [[db4ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip) :: nx,ny,nz,nq call me%destroy() nx = size(x,kind=ip) ny = size(y,kind=ip) nz = size(z,kind=ip) nq = size(q,kind=ip) call check_knot_vectors_sizes(nx=nx,kx=kx,tx=tx,& ny=ny,ky=ky,ty=ty,& nz=nz,kz=kz,tz=tz,& nq=nq,kq=kq,tq=tq,& iflag=iflag) if (iflag == 0_ip) then me%nx = nx me%ny = ny me%nz = nz me%nq = nq me%kx = kx me%ky = ky me%kz = kz me%kq = kq allocate(me%tx(nx+kx)) allocate(me%ty(ny+ky)) allocate(me%tz(nz+kz)) allocate(me%tq(nq+kq)) allocate(me%bcoef(nx,ny,nz,nq,1_ip)) allocate(me%work_val_1(ky,kz,kq)) allocate(me%work_val_2(kz,kq)) allocate(me%work_val_3(kq)) allocate(me%work_val_4(3_ip*max(kx,ky,kz,kq))) me%tx = tx me%ty = ty me%tz = tz me%tq = tq call db4ink(x,nx,y,ny,z,nz,q,nq,& fcn,& kx,ky,kz,kq,& 1_ip,& me%tx,me%ty,me%tz,me%tq,& me%bcoef(:,:,:,:,1),iflag) call me%set_extrap_flag(extrap) end if if (iflag==0_ip) me%nf = 1_ip me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_4d_specify_knots !***************************************************************************************** !***************************************************************************************** !> ! Evaluate a [[bspline_4d]] interpolate. This is a wrapper for [[db4val]]. pure subroutine evaluate_4d(me,xval,yval,zval,qval,idx,idy,idz,idq,f,iflag) implicit none class(bspline_4d),intent(inout) :: me real(wp),intent(in) :: xval !! \(x\) coordinate of evaluation point. real(wp),intent(in) :: yval !! \(y\) coordinate of evaluation point. real(wp),intent(in) :: zval !! \(z\) coordinate of evaluation point. real(wp),intent(in) :: qval !! \(q\) coordinate of evaluation point. integer(ip),intent(in) :: idx !! \(x\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idy !! \(y\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idz !! \(z\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idq !! \(q\) derivative of piecewise polynomial to evaluate. real(wp),intent(out) :: f !! interpolated value integer(ip),intent(out) :: iflag !! status flag (see [[db4val]]) if (me%initialized .and. me%nf/=1_ip) then f = 0.0_wp iflag = 6003_ip ! use the multiple-output version else if (me%initialized) then call db4val(xval,yval,zval,qval,& idx,idy,idz,idq,& me%tx,me%ty,me%tz,me%tq,& me%nx,me%ny,me%nz,me%nq,& me%kx,me%ky,me%kz,me%kq,& me%bcoef(:,:,:,:,1),f,iflag,& me%inbvx,me%inbvy,me%inbvz,me%inbvq,& me%iloy,me%iloz,me%iloq,& me%work_val_1,me%work_val_2,me%work_val_3,me%work_val_4,& extrap=me%extrap) else f = 0.0_wp iflag = 1_ip end if me%iflag = iflag end subroutine evaluate_4d !***************************************************************************************** !***************************************************************************************** !> ! It returns an empty [[bspline_5d]] type. Note that INITIALIZE still ! needs to be called before it can be used. ! Not really that useful except perhaps in some OpenMP applications. elemental function bspline_5d_constructor_empty() result(me) implicit none type(bspline_5d) :: me end function bspline_5d_constructor_empty !***************************************************************************************** !***************************************************************************************** !> ! Constructor for a [[bspline_5d]] type (auto knots). ! This is a wrapper for [[initialize_5d_auto_knots]]. pure function bspline_5d_constructor_auto_knots(x,y,z,q,r,fcn,kx,ky,kz,kq,kr,extrap) result(me) implicit none type(bspline_5d) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: r !! `(nr)` array of \(r\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr)` matrix of function values to interpolate. !! `fcn(i,j,k,l,m)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 2 \le k_q < n_q \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) !! ( \( 2 \le k_r < n_r \) ) !! (order = polynomial degree + 1) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) call initialize_5d_auto_knots(me,x,y,z,q,r,fcn,kx,ky,kz,kq,kr,me%iflag,extrap) end function bspline_5d_constructor_auto_knots !***************************************************************************************** !***************************************************************************************** !> ! Constructor for a [[bspline_5d]] type (user-specified knots). ! This is a wrapper for [[initialize_5d_specify_knots]]. pure function bspline_5d_constructor_specify_knots(x,y,z,q,r,fcn,& kx,ky,kz,kq,kr,& tx,ty,tz,tq,tr,extrap) result(me) implicit none type(bspline_5d) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: r !! `(nr)` array of \(r\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr)` matrix of function values to interpolate. !! `fcn(i,j,k,l,m)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 2 \le k_q < n_q \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) !! ( \( 2 \le k_r < n_r \) ) !! (order = polynomial degree + 1) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tr !! The `(nr+kr)` knots in the \(r\) direction !! for the spline interpolant. !! Must be non-decreasing. logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) call initialize_5d_specify_knots(me,x,y,z,q,r,fcn,kx,ky,kz,kq,kr,tx,ty,tz,tq,tr,me%iflag,extrap) end function bspline_5d_constructor_specify_knots !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_5d]] type (with automatically-computed knots). ! This is a wrapper for [[db5ink]]. pure subroutine initialize_5d_auto_knots(me,x,y,z,q,r,fcn,kx,ky,kz,kq,kr,iflag,extrap) implicit none class(bspline_5d),intent(inout) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: r !! `(nr)` array of \(r\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr)` matrix of function values to interpolate. !! `fcn(i,j,k,l,m)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 2 \le k_q < n_q \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) !! ( \( 2 \le k_r < n_r \) ) !! (order = polynomial degree + 1) integer(ip),intent(out) :: iflag !! status flag (see [[db5ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip) :: iknot integer(ip) :: nx,ny,nz,nq,nr call me%destroy() nx = size(x,kind=ip) ny = size(y,kind=ip) nz = size(z,kind=ip) nq = size(q,kind=ip) nr = size(r,kind=ip) me%nx = nx me%ny = ny me%nz = nz me%nq = nq me%nr = nr me%kx = kx me%ky = ky me%kz = kz me%kq = kq me%kr = kr allocate(me%tx(nx+kx)) allocate(me%ty(ny+ky)) allocate(me%tz(nz+kz)) allocate(me%tq(nq+kq)) allocate(me%tr(nr+kr)) allocate(me%bcoef(nx,ny,nz,nq,nr,1_ip)) allocate(me%work_val_1(ky,kz,kq,kr)) allocate(me%work_val_2(kz,kq,kr)) allocate(me%work_val_3(kq,kr)) allocate(me%work_val_4(kr)) allocate(me%work_val_5(3_ip*max(kx,ky,kz,kq,kr))) iknot = 0_ip !knot sequence chosen by db5ink call db5ink(x,nx,y,ny,z,nz,q,nq,r,nr,& fcn,& kx,ky,kz,kq,kr,& iknot,& me%tx,me%ty,me%tz,me%tq,me%tr,& me%bcoef(:,:,:,:,:,1),iflag) if (iflag==0_ip) then call me%set_extrap_flag(extrap) end if if (iflag==0_ip) me%nf = 1_ip me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_5d_auto_knots !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_5d]] type (with user-specified knots). ! This is a wrapper for [[db5ink]]. pure subroutine initialize_5d_specify_knots(me,x,y,z,q,r,fcn,& kx,ky,kz,kq,kr,& tx,ty,tz,tq,tr,iflag,extrap) implicit none class(bspline_5d),intent(inout) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: r !! `(nr)` array of \(r\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr)` matrix of function values to interpolate. !! `fcn(i,j,k,l,m)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 2 \le k_q < n_q \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) !! ( \( 2 \le k_r < n_r \) ) !! (order = polynomial degree + 1) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tr !! The `(nr+kr)` knots in the \(r\) direction !! for the spline interpolant. !! Must be non-decreasing. integer(ip),intent(out) :: iflag !! status flag (see [[db5ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip) :: nx,ny,nz,nq,nr call me%destroy() nx = size(x,kind=ip) ny = size(y,kind=ip) nz = size(z,kind=ip) nq = size(q,kind=ip) nr = size(r,kind=ip) call check_knot_vectors_sizes(nx=nx,kx=kx,tx=tx,& ny=ny,ky=ky,ty=ty,& nz=nz,kz=kz,tz=tz,& nq=nq,kq=kq,tq=tq,& nr=nr,kr=kr,tr=tr,& iflag=iflag) if (iflag == 0_ip) then me%nx = nx me%ny = ny me%nz = nz me%nq = nq me%nr = nr me%kx = kx me%ky = ky me%kz = kz me%kq = kq me%kr = kr allocate(me%tx(nx+kx)) allocate(me%ty(ny+ky)) allocate(me%tz(nz+kz)) allocate(me%tq(nq+kq)) allocate(me%tr(nr+kr)) allocate(me%bcoef(nx,ny,nz,nq,nr,1_ip)) allocate(me%work_val_1(ky,kz,kq,kr)) allocate(me%work_val_2(kz,kq,kr)) allocate(me%work_val_3(kq,kr)) allocate(me%work_val_4(kr)) allocate(me%work_val_5(3_ip*max(kx,ky,kz,kq,kr))) me%tx = tx me%ty = ty me%tz = tz me%tq = tq me%tr = tr call db5ink(x,nx,y,ny,z,nz,q,nq,r,nr,& fcn,& kx,ky,kz,kq,kr,& 1_ip,& me%tx,me%ty,me%tz,me%tq,me%tr,& me%bcoef(:,:,:,:,:,1),iflag) call me%set_extrap_flag(extrap) end if if (iflag==0_ip) me%nf = 1_ip me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_5d_specify_knots !***************************************************************************************** !***************************************************************************************** !> ! Evaluate a [[bspline_5d]] interpolate. This is a wrapper for [[db5val]]. pure subroutine evaluate_5d(me,xval,yval,zval,qval,rval,idx,idy,idz,idq,idr,f,iflag) implicit none class(bspline_5d),intent(inout) :: me real(wp),intent(in) :: xval !! \(x\) coordinate of evaluation point. real(wp),intent(in) :: yval !! \(y\) coordinate of evaluation point. real(wp),intent(in) :: zval !! \(z\) coordinate of evaluation point. real(wp),intent(in) :: qval !! \(q\) coordinate of evaluation point. real(wp),intent(in) :: rval !! \(r\) coordinate of evaluation point. integer(ip),intent(in) :: idx !! \(x\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idy !! \(y\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idz !! \(z\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idq !! \(q\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idr !! \(r\) derivative of piecewise polynomial to evaluate. real(wp),intent(out) :: f !! interpolated value integer(ip),intent(out) :: iflag !! status flag (see [[db5val]]) if (me%initialized .and. me%nf/=1_ip) then f = 0.0_wp iflag = 6003_ip ! use the multiple-output version else if (me%initialized) then call db5val(xval,yval,zval,qval,rval,& idx,idy,idz,idq,idr,& me%tx,me%ty,me%tz,me%tq,me%tr,& me%nx,me%ny,me%nz,me%nq,me%nr,& me%kx,me%ky,me%kz,me%kq,me%kr,& me%bcoef(:,:,:,:,:,1),f,iflag,& me%inbvx,me%inbvy,me%inbvz,me%inbvq,me%inbvr,& me%iloy,me%iloz,me%iloq,me%ilor,& me%work_val_1,me%work_val_2,me%work_val_3,me%work_val_4,me%work_val_5,& extrap=me%extrap) else f = 0.0_wp iflag = 1_ip end if me%iflag = iflag end subroutine evaluate_5d !***************************************************************************************** !***************************************************************************************** !> ! It returns an empty [[bspline_6d]] type. Note that INITIALIZE still ! needs to be called before it can be used. ! Not really that useful except perhaps in some OpenMP applications. elemental function bspline_6d_constructor_empty() result(me) implicit none type(bspline_6d) :: me end function bspline_6d_constructor_empty !***************************************************************************************** !***************************************************************************************** !> ! Constructor for a [[bspline_6d]] type (auto knots). ! This is a wrapper for [[initialize_6d_auto_knots]]. pure function bspline_6d_constructor_auto_knots(x,y,z,q,r,s,fcn,& kx,ky,kz,kq,kr,ks,extrap) result(me) implicit none type(bspline_6d) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: r !! `(nr)` array of \(r\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: s !! `(ns)` array of \(s\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr,ns)` matrix of function values to interpolate. !! `fcn(i,j,k,l,m,n)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`,`s(n)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 2 \le k_q < n_q \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) !! ( \( 2 \le k_r < n_r \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ks !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) call initialize_6d_auto_knots(me,x,y,z,q,r,s,fcn,kx,ky,kz,kq,kr,ks,me%iflag,extrap) end function bspline_6d_constructor_auto_knots !***************************************************************************************** !***************************************************************************************** !> ! Constructor for a [[bspline_6d]] type (user-specified knots). ! This is a wrapper for [[initialize_6d_specify_knots]]. pure function bspline_6d_constructor_specify_knots(x,y,z,q,r,s,fcn,& kx,ky,kz,kq,kr,ks,& tx,ty,tz,tq,tr,ts,extrap) result(me) implicit none type(bspline_6d) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: r !! `(nr)` array of \(r\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: s !! `(ns)` array of \(s\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr,ns)` matrix of function values to interpolate. !! `fcn(i,j,k,l,m,n)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`,`s(n)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 2 \le k_q < n_q \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) !! ( \( 2 \le k_r < n_r \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ks !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tr !! The `(nr+kr)` knots in the \(r\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ts !! The `(ns+ks)` knots in the \(s\) direction !! for the spline interpolant. !! Must be non-decreasing. logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) call initialize_6d_specify_knots(me,x,y,z,q,r,s,fcn,& kx,ky,kz,kq,kr,ks,& tx,ty,tz,tq,tr,ts,me%iflag,extrap) end function bspline_6d_constructor_specify_knots !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_6d]] type (with automatically-computed knots). ! This is a wrapper for [[db6ink]]. pure subroutine initialize_6d_auto_knots(me,x,y,z,q,r,s,fcn,& kx,ky,kz,kq,kr,ks,iflag,extrap) implicit none class(bspline_6d),intent(inout) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: r !! `(nr)` array of \(r\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: s !! `(ns)` array of \(s\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr,ns)` matrix of function values to interpolate. !! `fcn(i,j,k,l,m,n)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`,`s(n)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 2 \le k_q < n_q \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) !! ( \( 2 \le k_r < n_r \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ks !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(out) :: iflag !! status flag (see [[db6ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip) :: iknot integer(ip) :: nx,ny,nz,nq,nr,ns call me%destroy() nx = size(x,kind=ip) ny = size(y,kind=ip) nz = size(z,kind=ip) nq = size(q,kind=ip) nr = size(r,kind=ip) ns = size(s,kind=ip) me%nx = nx me%ny = ny me%nz = nz me%nq = nq me%nr = nr me%ns = ns me%kx = kx me%ky = ky me%kz = kz me%kq = kq me%kr = kr me%ks = ks allocate(me%tx(nx+kx)) allocate(me%ty(ny+ky)) allocate(me%tz(nz+kz)) allocate(me%tq(nq+kq)) allocate(me%tr(nr+kr)) allocate(me%ts(ns+ks)) allocate(me%bcoef(nx,ny,nz,nq,nr,ns,1_ip)) allocate(me%work_val_1(ky,kz,kq,kr,ks)) allocate(me%work_val_2(kz,kq,kr,ks)) allocate(me%work_val_3(kq,kr,ks)) allocate(me%work_val_4(kr,ks)) allocate(me%work_val_5(ks)) allocate(me%work_val_6(3_ip*max(kx,ky,kz,kq,kr,ks))) iknot = 0_ip !knot sequence chosen by db6ink call db6ink(x,nx,y,ny,z,nz,q,nq,r,nr,s,ns,& fcn,& kx,ky,kz,kq,kr,ks,& iknot,& me%tx,me%ty,me%tz,me%tq,me%tr,me%ts,& me%bcoef(:,:,:,:,:,:,1),iflag) if (iflag==0_ip) then call me%set_extrap_flag(extrap) end if if (iflag==0_ip) me%nf = 1_ip me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_6d_auto_knots !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_6d]] type (with user-specified knots). ! This is a wrapper for [[db6ink]]. pure subroutine initialize_6d_specify_knots(me,x,y,z,q,r,s,fcn,& kx,ky,kz,kq,kr,ks,& tx,ty,tz,tq,tr,ts,iflag,extrap) implicit none class(bspline_6d),intent(inout) :: me real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: r !! `(nr)` array of \(r\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: s !! `(ns)` array of \(s\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr,ns)` matrix of function values to interpolate. !! `fcn(i,j,k,l,m,n)` should contain the function value at the !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`,`s(n)`) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 2 \le k_x < n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 2 \le k_y < n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 2 \le k_q < n_q \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) !! ( \( 2 \le k_r < n_r \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ks !! The order of spline pieces in \(z\) !! ( \( 2 \le k_z < n_z \) ) !! (order = polynomial degree + 1) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tr !! The `(nr+kr)` knots in the \(r\) direction !! for the spline interpolant. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ts !! The `(ns+ks)` knots in the \(s\) direction !! for the spline interpolant. !! Must be non-decreasing. integer(ip),intent(out) :: iflag !! status flag (see [[db6ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip) :: nx,ny,nz,nq,nr,ns call me%destroy() nx = size(x,kind=ip) ny = size(y,kind=ip) nz = size(z,kind=ip) nq = size(q,kind=ip) nr = size(r,kind=ip) ns = size(s,kind=ip) call check_knot_vectors_sizes(nx=nx,kx=kx,tx=tx,& ny=ny,ky=ky,ty=ty,& nz=nz,kz=kz,tz=tz,& nq=nq,kq=kq,tq=tq,& nr=nr,kr=kr,tr=tr,& ns=ns,ks=ks,ts=ts,& iflag=iflag) if (iflag == 0_ip) then me%nx = nx me%ny = ny me%nz = nz me%nq = nq me%nr = nr me%ns = ns me%kx = kx me%ky = ky me%kz = kz me%kq = kq me%kr = kr me%ks = ks allocate(me%tx(nx+kx)) allocate(me%ty(ny+ky)) allocate(me%tz(nz+kz)) allocate(me%tq(nq+kq)) allocate(me%tr(nr+kr)) allocate(me%ts(ns+ks)) allocate(me%bcoef(nx,ny,nz,nq,nr,ns,1_ip)) allocate(me%work_val_1(ky,kz,kq,kr,ks)) allocate(me%work_val_2(kz,kq,kr,ks)) allocate(me%work_val_3(kq,kr,ks)) allocate(me%work_val_4(kr,ks)) allocate(me%work_val_5(ks)) allocate(me%work_val_6(3_ip*max(kx,ky,kz,kq,kr,ks))) me%tx = tx me%ty = ty me%tz = tz me%tq = tq me%tr = tr me%ts = ts call db6ink(x,nx,y,ny,z,nz,q,nq,r,nr,s,ns,& fcn,& kx,ky,kz,kq,kr,ks,& 1_ip,& me%tx,me%ty,me%tz,me%tq,me%tr,me%ts,& me%bcoef(:,:,:,:,:,:,1),iflag) call me%set_extrap_flag(extrap) end if if (iflag==0_ip) me%nf = 1_ip me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_6d_specify_knots !***************************************************************************************** !***************************************************************************************** !> ! Evaluate a [[bspline_6d]] interpolate. This is a wrapper for [[db6val]]. pure subroutine evaluate_6d(me,xval,yval,zval,qval,rval,sval,idx,idy,idz,idq,idr,ids,f,iflag) implicit none class(bspline_6d),intent(inout) :: me real(wp),intent(in) :: xval !! \(x\) coordinate of evaluation point. real(wp),intent(in) :: yval !! \(y\) coordinate of evaluation point. real(wp),intent(in) :: zval !! \(z\) coordinate of evaluation point. real(wp),intent(in) :: qval !! \(q\) coordinate of evaluation point. real(wp),intent(in) :: rval !! \(r\) coordinate of evaluation point. real(wp),intent(in) :: sval !! \(s\) coordinate of evaluation point. integer(ip),intent(in) :: idx !! \(x\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idy !! \(y\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idz !! \(z\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idq !! \(q\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idr !! \(r\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: ids !! \(s\) derivative of piecewise polynomial to evaluate. real(wp),intent(out) :: f !! interpolated value integer(ip),intent(out) :: iflag !! status flag (see [[db6val]]) if (me%initialized .and. me%nf/=1_ip) then f = 0.0_wp iflag = 6003_ip ! use the multiple-output version else if (me%initialized) then call db6val(xval,yval,zval,qval,rval,sval,& idx,idy,idz,idq,idr,ids,& me%tx,me%ty,me%tz,me%tq,me%tr,me%ts,& me%nx,me%ny,me%nz,me%nq,me%nr,me%ns,& me%kx,me%ky,me%kz,me%kq,me%kr,me%ks,& me%bcoef(:,:,:,:,:,:,1),f,iflag,& me%inbvx,me%inbvy,me%inbvz,me%inbvq,me%inbvr,me%inbvs,& me%iloy,me%iloz,me%iloq,me%ilor,me%ilos,& me%work_val_1,me%work_val_2,me%work_val_3,me%work_val_4,me%work_val_5,me%work_val_6,& extrap=me%extrap) else f = 0.0_wp iflag = 1_ip end if me%iflag = iflag end subroutine evaluate_6d !***************************************************************************************** !***************************************************************************************** !> ! Error checks for the user-specified knot vector sizes. ! !@note If more than one is the wrong size, then the `iflag` error code will ! correspond to the one with the highest rank. pure subroutine check_knot_vectors_sizes(nx,ny,nz,nq,nr,ns,& kx,ky,kz,kq,kr,ks,& tx,ty,tz,tq,tr,ts,iflag) implicit none integer(ip),intent(in),optional :: nx integer(ip),intent(in),optional :: ny integer(ip),intent(in),optional :: nz integer(ip),intent(in),optional :: nq integer(ip),intent(in),optional :: nr integer(ip),intent(in),optional :: ns integer(ip),intent(in),optional :: kx integer(ip),intent(in),optional :: ky integer(ip),intent(in),optional :: kz integer(ip),intent(in),optional :: kq integer(ip),intent(in),optional :: kr integer(ip),intent(in),optional :: ks real(wp),dimension(:),intent(in),optional :: tx real(wp),dimension(:),intent(in),optional :: ty real(wp),dimension(:),intent(in),optional :: tz real(wp),dimension(:),intent(in),optional :: tq real(wp),dimension(:),intent(in),optional :: tr real(wp),dimension(:),intent(in),optional :: ts integer(ip),intent(out) :: iflag !! 0 if everything is OK iflag = 0_ip if (present(nx) .and. present(kx) .and. present(tx)) then if (size(tx,kind=ip)/=(nx+kx)) then iflag = 501_ip ! tx is not the correct size (nx+kx) end if end if if (present(ny) .and. present(ky) .and. present(ty)) then if (size(ty,kind=ip)/=(ny+ky)) then iflag = 502_ip ! ty is not the correct size (ny+ky) end if end if if (present(nz) .and. present(kz) .and. present(tz)) then if (size(tz,kind=ip)/=(nz+kz)) then iflag = 503_ip ! tz is not the correct size (nz+kz) end if end if if (present(nq) .and. present(kq) .and. present(tq)) then if (size(tq,kind=ip)/=(nq+kq)) then iflag = 504_ip ! tq is not the correct size (nq+kq) end if end if if (present(nr) .and. present(kr) .and. present(tr)) then if (size(tr,kind=ip)/=(nr+kr)) then iflag = 505_ip ! tr is not the correct size (nr+kr) end if end if if (present(ns) .and. present(ks) .and. present(ts)) then if (size(ts,kind=ip)/=(ns+ks)) then iflag = 506_ip ! ts is not the correct size (ns+ks) end if end if end subroutine check_knot_vectors_sizes !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_1d]] type directly from its knots and b-spline ! coefficients (for example, those saved from another object with ! `get_knots` and `get_coefficients`), without interpolating any data. ! ! The number of coefficients in each direction is taken from the ! shape of `bcoef`. pure subroutine initialize_from_coefficients_1d(me,tx,bcoef,kx,iflag,extrap) implicit none class(bspline_1d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: bcoef !! `(nx)` array of b-spline coefficients integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 1 \le k_x \le n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(out) :: iflag !! status flag: !! !! * 0: no errors !! * 5001-5005: invalid knots or coefficients !! (see [[get_status_message]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(1) :: n !! number of coefficients in each direction call me%destroy() n = shape(bcoef,kind=ip) call check_coefficients(n(1),kx,tx,iflag) if (iflag == 0_ip) then me%nx = n(1) me%kx = kx allocate(me%bcoef(n(1),1_ip)) ! (assigned below: no temporary array) me%bcoef(:,1) = bcoef me%nf = 1_ip me%tx = tx allocate(me%work_val_1(3_ip*kx)) call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_from_coefficients_1d !***************************************************************************************** !***************************************************************************************** !> ! Returns the knot vectors of a [[bspline_1d]]. pure subroutine get_knots_1d(me,tx,iflag) implicit none class(bspline_1d),intent(inout) :: me !! the spline real(wp),dimension(:),allocatable,intent(out) :: tx !! the `(nx+kx)` knots in the \(x\) direction integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized) then tx = me%tx iflag = 0_ip else allocate(tx(0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_knots_1d !***************************************************************************************** !***************************************************************************************** !> ! Returns the b-spline coefficients of a [[bspline_1d]]. ! These are the values of the control points (see [[greville_abscissae]] ! for their locations). pure subroutine get_coefficients_1d(me,bcoef,iflag) implicit none class(bspline_1d),intent(inout) :: me !! the spline real(wp),dimension(:),allocatable,intent(out) :: bcoef !! `(nx)` array of b-spline coefficients integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized .and. me%nf/=1_ip) then iflag = 6003_ip ! use the multiple-output version allocate(bcoef(0)) else if (me%initialized) then bcoef = me%bcoef(:,1) iflag = 0_ip else allocate(bcoef(0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_coefficients_1d !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_2d]] type directly from its knots and b-spline ! coefficients (for example, those saved from another object with ! `get_knots` and `get_coefficients`), without interpolating any data. ! ! The number of coefficients in each direction is taken from the ! shape of `bcoef`. pure subroutine initialize_from_coefficients_2d(me,tx,ty,bcoef,kx,ky,iflag,extrap) implicit none class(bspline_2d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. !! Must be non-decreasing. real(wp),dimension(:,:),intent(in) :: bcoef !! `(nx,ny)` array of b-spline coefficients integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 1 \le k_x \le n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 1 \le k_y \le n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(out) :: iflag !! status flag: !! !! * 0: no errors !! * 5001-5005: invalid knots or coefficients !! (see [[get_status_message]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(2) :: n !! number of coefficients in each direction call me%destroy() n = shape(bcoef,kind=ip) call check_coefficients(n(1),kx,tx,iflag) if (iflag==0_ip) call check_coefficients(n(2),ky,ty,iflag) if (iflag == 0_ip) then me%nx = n(1) me%ny = n(2) me%kx = kx me%ky = ky allocate(me%bcoef(n(1),n(2),1_ip)) ! (assigned below: no temporary array) me%bcoef(:,:,1) = bcoef me%nf = 1_ip me%tx = tx me%ty = ty allocate(me%work_val_1(ky)) allocate(me%work_val_2(3_ip*max(kx,ky))) call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_from_coefficients_2d !***************************************************************************************** !***************************************************************************************** !> ! Returns the knot vectors of a [[bspline_2d]]. pure subroutine get_knots_2d(me,tx,ty,iflag) implicit none class(bspline_2d),intent(inout) :: me !! the spline real(wp),dimension(:),allocatable,intent(out) :: tx !! the `(nx+kx)` knots in the \(x\) direction real(wp),dimension(:),allocatable,intent(out) :: ty !! the `(ny+ky)` knots in the \(y\) direction integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized) then tx = me%tx ty = me%ty iflag = 0_ip else allocate(tx(0)) allocate(ty(0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_knots_2d !***************************************************************************************** !***************************************************************************************** !> ! Returns the b-spline coefficients of a [[bspline_2d]]. ! These are the values of the control points (see [[greville_abscissae]] ! for their locations). pure subroutine get_coefficients_2d(me,bcoef,iflag) implicit none class(bspline_2d),intent(inout) :: me !! the spline real(wp),dimension(:,:),allocatable,intent(out) :: bcoef !! `(nx,ny)` array of b-spline coefficients integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized .and. me%nf/=1_ip) then iflag = 6003_ip ! use the multiple-output version allocate(bcoef(0,0)) else if (me%initialized) then bcoef = me%bcoef(:,:,1) iflag = 0_ip else allocate(bcoef(0,0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_coefficients_2d !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_3d]] type directly from its knots and b-spline ! coefficients (for example, those saved from another object with ! `get_knots` and `get_coefficients`), without interpolating any data. ! ! The number of coefficients in each direction is taken from the ! shape of `bcoef`. pure subroutine initialize_from_coefficients_3d(me,tx,ty,tz,bcoef,kx,ky,kz,iflag,extrap) implicit none class(bspline_3d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction. !! Must be non-decreasing. real(wp),dimension(:,:,:),intent(in) :: bcoef !! `(nx,ny,nz)` array of b-spline coefficients integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 1 \le k_x \le n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 1 \le k_y \le n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 1 \le k_z \le n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(out) :: iflag !! status flag: !! !! * 0: no errors !! * 5001-5005: invalid knots or coefficients !! (see [[get_status_message]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(3) :: n !! number of coefficients in each direction call me%destroy() n = shape(bcoef,kind=ip) call check_coefficients(n(1),kx,tx,iflag) if (iflag==0_ip) call check_coefficients(n(2),ky,ty,iflag) if (iflag==0_ip) call check_coefficients(n(3),kz,tz,iflag) if (iflag == 0_ip) then me%nx = n(1) me%ny = n(2) me%nz = n(3) me%kx = kx me%ky = ky me%kz = kz allocate(me%bcoef(n(1),n(2),n(3),1_ip)) ! (assigned below: no temporary array) me%bcoef(:,:,:,1) = bcoef me%nf = 1_ip me%tx = tx me%ty = ty me%tz = tz allocate(me%work_val_1(ky,kz)) allocate(me%work_val_2(kz)) allocate(me%work_val_3(3_ip*max(kx,ky,kz))) call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_from_coefficients_3d !***************************************************************************************** !***************************************************************************************** !> ! Returns the knot vectors of a [[bspline_3d]]. pure subroutine get_knots_3d(me,tx,ty,tz,iflag) implicit none class(bspline_3d),intent(inout) :: me !! the spline real(wp),dimension(:),allocatable,intent(out) :: tx !! the `(nx+kx)` knots in the \(x\) direction real(wp),dimension(:),allocatable,intent(out) :: ty !! the `(ny+ky)` knots in the \(y\) direction real(wp),dimension(:),allocatable,intent(out) :: tz !! the `(nz+kz)` knots in the \(z\) direction integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized) then tx = me%tx ty = me%ty tz = me%tz iflag = 0_ip else allocate(tx(0)) allocate(ty(0)) allocate(tz(0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_knots_3d !***************************************************************************************** !***************************************************************************************** !> ! Returns the b-spline coefficients of a [[bspline_3d]]. ! These are the values of the control points (see [[greville_abscissae]] ! for their locations). pure subroutine get_coefficients_3d(me,bcoef,iflag) implicit none class(bspline_3d),intent(inout) :: me !! the spline real(wp),dimension(:,:,:),allocatable,intent(out) :: bcoef !! `(nx,ny,nz)` array of b-spline coefficients integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized .and. me%nf/=1_ip) then iflag = 6003_ip ! use the multiple-output version allocate(bcoef(0,0,0)) else if (me%initialized) then bcoef = me%bcoef(:,:,:,1) iflag = 0_ip else allocate(bcoef(0,0,0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_coefficients_3d !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_4d]] type directly from its knots and b-spline ! coefficients (for example, those saved from another object with ! `get_knots` and `get_coefficients`), without interpolating any data. ! ! The number of coefficients in each direction is taken from the ! shape of `bcoef`. pure subroutine initialize_from_coefficients_4d(me,tx,ty,tz,tq,bcoef,kx,ky,kz,kq,iflag,extrap) implicit none class(bspline_4d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction. !! Must be non-decreasing. real(wp),dimension(:,:,:,:),intent(in) :: bcoef !! `(nx,ny,nz,nq)` array of b-spline coefficients integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 1 \le k_x \le n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 1 \le k_y \le n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 1 \le k_z \le n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 1 \le k_q \le n_q \) ) !! (order = polynomial degree + 1) integer(ip),intent(out) :: iflag !! status flag: !! !! * 0: no errors !! * 5001-5005: invalid knots or coefficients !! (see [[get_status_message]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(4) :: n !! number of coefficients in each direction call me%destroy() n = shape(bcoef,kind=ip) call check_coefficients(n(1),kx,tx,iflag) if (iflag==0_ip) call check_coefficients(n(2),ky,ty,iflag) if (iflag==0_ip) call check_coefficients(n(3),kz,tz,iflag) if (iflag==0_ip) call check_coefficients(n(4),kq,tq,iflag) if (iflag == 0_ip) then me%nx = n(1) me%ny = n(2) me%nz = n(3) me%nq = n(4) me%kx = kx me%ky = ky me%kz = kz me%kq = kq allocate(me%bcoef(n(1),n(2),n(3),n(4),1_ip)) ! (assigned below: no temporary array) me%bcoef(:,:,:,:,1) = bcoef me%nf = 1_ip me%tx = tx me%ty = ty me%tz = tz me%tq = tq allocate(me%work_val_1(ky,kz,kq)) allocate(me%work_val_2(kz,kq)) allocate(me%work_val_3(kq)) allocate(me%work_val_4(3_ip*max(kx,ky,kz,kq))) call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_from_coefficients_4d !***************************************************************************************** !***************************************************************************************** !> ! Returns the knot vectors of a [[bspline_4d]]. pure subroutine get_knots_4d(me,tx,ty,tz,tq,iflag) implicit none class(bspline_4d),intent(inout) :: me !! the spline real(wp),dimension(:),allocatable,intent(out) :: tx !! the `(nx+kx)` knots in the \(x\) direction real(wp),dimension(:),allocatable,intent(out) :: ty !! the `(ny+ky)` knots in the \(y\) direction real(wp),dimension(:),allocatable,intent(out) :: tz !! the `(nz+kz)` knots in the \(z\) direction real(wp),dimension(:),allocatable,intent(out) :: tq !! the `(nq+kq)` knots in the \(q\) direction integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized) then tx = me%tx ty = me%ty tz = me%tz tq = me%tq iflag = 0_ip else allocate(tx(0)) allocate(ty(0)) allocate(tz(0)) allocate(tq(0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_knots_4d !***************************************************************************************** !***************************************************************************************** !> ! Returns the b-spline coefficients of a [[bspline_4d]]. ! These are the values of the control points (see [[greville_abscissae]] ! for their locations). pure subroutine get_coefficients_4d(me,bcoef,iflag) implicit none class(bspline_4d),intent(inout) :: me !! the spline real(wp),dimension(:,:,:,:),allocatable,intent(out) :: bcoef !! `(nx,ny,nz,nq)` array of b-spline coefficients integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized .and. me%nf/=1_ip) then iflag = 6003_ip ! use the multiple-output version allocate(bcoef(0,0,0,0)) else if (me%initialized) then bcoef = me%bcoef(:,:,:,:,1) iflag = 0_ip else allocate(bcoef(0,0,0,0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_coefficients_4d !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_5d]] type directly from its knots and b-spline ! coefficients (for example, those saved from another object with ! `get_knots` and `get_coefficients`), without interpolating any data. ! ! The number of coefficients in each direction is taken from the ! shape of `bcoef`. pure subroutine initialize_from_coefficients_5d(me,tx,ty,tz,tq,tr,bcoef,kx,ky,kz,kq,kr,iflag,extrap) implicit none class(bspline_5d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tr !! The `(nr+kr)` knots in the \(r\) direction. !! Must be non-decreasing. real(wp),dimension(:,:,:,:,:),intent(in) :: bcoef !! `(nx,ny,nz,nq,nr)` array of b-spline coefficients integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 1 \le k_x \le n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 1 \le k_y \le n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 1 \le k_z \le n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 1 \le k_q \le n_q \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) !! ( \( 1 \le k_r \le n_r \) ) !! (order = polynomial degree + 1) integer(ip),intent(out) :: iflag !! status flag: !! !! * 0: no errors !! * 5001-5005: invalid knots or coefficients !! (see [[get_status_message]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(5) :: n !! number of coefficients in each direction call me%destroy() n = shape(bcoef,kind=ip) call check_coefficients(n(1),kx,tx,iflag) if (iflag==0_ip) call check_coefficients(n(2),ky,ty,iflag) if (iflag==0_ip) call check_coefficients(n(3),kz,tz,iflag) if (iflag==0_ip) call check_coefficients(n(4),kq,tq,iflag) if (iflag==0_ip) call check_coefficients(n(5),kr,tr,iflag) if (iflag == 0_ip) then me%nx = n(1) me%ny = n(2) me%nz = n(3) me%nq = n(4) me%nr = n(5) me%kx = kx me%ky = ky me%kz = kz me%kq = kq me%kr = kr allocate(me%bcoef(n(1),n(2),n(3),n(4),n(5),1_ip)) ! (assigned below: no temporary array) me%bcoef(:,:,:,:,:,1) = bcoef me%nf = 1_ip me%tx = tx me%ty = ty me%tz = tz me%tq = tq me%tr = tr allocate(me%work_val_1(ky,kz,kq,kr)) allocate(me%work_val_2(kz,kq,kr)) allocate(me%work_val_3(kq,kr)) allocate(me%work_val_4(kr)) allocate(me%work_val_5(3_ip*max(kx,ky,kz,kq,kr))) call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_from_coefficients_5d !***************************************************************************************** !***************************************************************************************** !> ! Returns the knot vectors of a [[bspline_5d]]. pure subroutine get_knots_5d(me,tx,ty,tz,tq,tr,iflag) implicit none class(bspline_5d),intent(inout) :: me !! the spline real(wp),dimension(:),allocatable,intent(out) :: tx !! the `(nx+kx)` knots in the \(x\) direction real(wp),dimension(:),allocatable,intent(out) :: ty !! the `(ny+ky)` knots in the \(y\) direction real(wp),dimension(:),allocatable,intent(out) :: tz !! the `(nz+kz)` knots in the \(z\) direction real(wp),dimension(:),allocatable,intent(out) :: tq !! the `(nq+kq)` knots in the \(q\) direction real(wp),dimension(:),allocatable,intent(out) :: tr !! the `(nr+kr)` knots in the \(r\) direction integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized) then tx = me%tx ty = me%ty tz = me%tz tq = me%tq tr = me%tr iflag = 0_ip else allocate(tx(0)) allocate(ty(0)) allocate(tz(0)) allocate(tq(0)) allocate(tr(0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_knots_5d !***************************************************************************************** !***************************************************************************************** !> ! Returns the b-spline coefficients of a [[bspline_5d]]. ! These are the values of the control points (see [[greville_abscissae]] ! for their locations). pure subroutine get_coefficients_5d(me,bcoef,iflag) implicit none class(bspline_5d),intent(inout) :: me !! the spline real(wp),dimension(:,:,:,:,:),allocatable,intent(out) :: bcoef !! `(nx,ny,nz,nq,nr)` array of b-spline coefficients integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized .and. me%nf/=1_ip) then iflag = 6003_ip ! use the multiple-output version allocate(bcoef(0,0,0,0,0)) else if (me%initialized) then bcoef = me%bcoef(:,:,:,:,:,1) iflag = 0_ip else allocate(bcoef(0,0,0,0,0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_coefficients_5d !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_6d]] type directly from its knots and b-spline ! coefficients (for example, those saved from another object with ! `get_knots` and `get_coefficients`), without interpolating any data. ! ! The number of coefficients in each direction is taken from the ! shape of `bcoef`. pure subroutine initialize_from_coefficients_6d(me,tx,ty,tz,tq,tr,ts,bcoef,kx,ky,kz,kq,kr,ks,iflag,extrap) implicit none class(bspline_6d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: tr !! The `(nr+kr)` knots in the \(r\) direction. !! Must be non-decreasing. real(wp),dimension(:),intent(in) :: ts !! The `(ns+ks)` knots in the \(s\) direction. !! Must be non-decreasing. real(wp),dimension(:,:,:,:,:,:),intent(in) :: bcoef !! `(nx,ny,nz,nq,nr,ns)` array of b-spline coefficients integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) !! ( \( 1 \le k_x \le n_x \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) !! ( \( 1 \le k_y \le n_y \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) !! ( \( 1 \le k_z \le n_z \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) !! ( \( 1 \le k_q \le n_q \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) !! ( \( 1 \le k_r \le n_r \) ) !! (order = polynomial degree + 1) integer(ip),intent(in) :: ks !! The order of spline pieces in \(s\) !! ( \( 1 \le k_s \le n_s \) ) !! (order = polynomial degree + 1) integer(ip),intent(out) :: iflag !! status flag: !! !! * 0: no errors !! * 5001-5005: invalid knots or coefficients !! (see [[get_status_message]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(6) :: n !! number of coefficients in each direction call me%destroy() n = shape(bcoef,kind=ip) call check_coefficients(n(1),kx,tx,iflag) if (iflag==0_ip) call check_coefficients(n(2),ky,ty,iflag) if (iflag==0_ip) call check_coefficients(n(3),kz,tz,iflag) if (iflag==0_ip) call check_coefficients(n(4),kq,tq,iflag) if (iflag==0_ip) call check_coefficients(n(5),kr,tr,iflag) if (iflag==0_ip) call check_coefficients(n(6),ks,ts,iflag) if (iflag == 0_ip) then me%nx = n(1) me%ny = n(2) me%nz = n(3) me%nq = n(4) me%nr = n(5) me%ns = n(6) me%kx = kx me%ky = ky me%kz = kz me%kq = kq me%kr = kr me%ks = ks allocate(me%bcoef(n(1),n(2),n(3),n(4),n(5),n(6),1_ip)) ! (assigned below: no temporary array) me%bcoef(:,:,:,:,:,:,1) = bcoef me%nf = 1_ip me%tx = tx me%ty = ty me%tz = tz me%tq = tq me%tr = tr me%ts = ts allocate(me%work_val_1(ky,kz,kq,kr,ks)) allocate(me%work_val_2(kz,kq,kr,ks)) allocate(me%work_val_3(kq,kr,ks)) allocate(me%work_val_4(kr,ks)) allocate(me%work_val_5(ks)) allocate(me%work_val_6(3_ip*max(kx,ky,kz,kq,kr,ks))) call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_from_coefficients_6d !***************************************************************************************** !***************************************************************************************** !> ! Returns the knot vectors of a [[bspline_6d]]. pure subroutine get_knots_6d(me,tx,ty,tz,tq,tr,ts,iflag) implicit none class(bspline_6d),intent(inout) :: me !! the spline real(wp),dimension(:),allocatable,intent(out) :: tx !! the `(nx+kx)` knots in the \(x\) direction real(wp),dimension(:),allocatable,intent(out) :: ty !! the `(ny+ky)` knots in the \(y\) direction real(wp),dimension(:),allocatable,intent(out) :: tz !! the `(nz+kz)` knots in the \(z\) direction real(wp),dimension(:),allocatable,intent(out) :: tq !! the `(nq+kq)` knots in the \(q\) direction real(wp),dimension(:),allocatable,intent(out) :: tr !! the `(nr+kr)` knots in the \(r\) direction real(wp),dimension(:),allocatable,intent(out) :: ts !! the `(ns+ks)` knots in the \(s\) direction integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized) then tx = me%tx ty = me%ty tz = me%tz tq = me%tq tr = me%tr ts = me%ts iflag = 0_ip else allocate(tx(0)) allocate(ty(0)) allocate(tz(0)) allocate(tq(0)) allocate(tr(0)) allocate(ts(0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_knots_6d !***************************************************************************************** !***************************************************************************************** !> ! Returns the b-spline coefficients of a [[bspline_6d]]. ! These are the values of the control points (see [[greville_abscissae]] ! for their locations). pure subroutine get_coefficients_6d(me,bcoef,iflag) implicit none class(bspline_6d),intent(inout) :: me !! the spline real(wp),dimension(:,:,:,:,:,:),allocatable,intent(out) :: bcoef !! `(nx,ny,nz,nq,nr,ns)` array of b-spline coefficients integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized .and. me%nf/=1_ip) then iflag = 6003_ip ! use the multiple-output version allocate(bcoef(0,0,0,0,0,0)) else if (me%initialized) then bcoef = me%bcoef(:,:,:,:,:,:,1) iflag = 0_ip else allocate(bcoef(0,0,0,0,0,0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_coefficients_6d !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_1d]] type with multiple outputs: `nf` functions of ! the same inputs, on the same grid. ! The last dimension of `fcn` is the output index. ! This is a wrapper for [[db1ink]] (multiple-output version). pure subroutine initialize_1d_auto_knots_multi(me,x,fcn,kx,iflag,extrap) implicit none class(bspline_1d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:,:),intent(in) :: fcn !! `(nx,nf)` array of function values to interpolate !! (the last dimension is the output index) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(out) :: iflag !! status flag (see [[db1ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(1) :: n !! number of points in each direction integer(ip) :: nf !! number of outputs call me%destroy() n(1) = size(x,kind=ip) nf = size(fcn,2,kind=ip) me%nx = n(1) me%kx = kx allocate(me%tx(n(1)+kx)) allocate(me%bcoef(n(1),nf)) allocate(me%work_val_1(3_ip*kx)) call db1ink(x,n(1),fcn,kx,0_ip,me%tx,me%bcoef,iflag) if (iflag==0_ip) then me%nf = nf call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_1d_auto_knots_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_1d]] type with multiple outputs: `nf` functions of ! the same inputs, on the same grid (with user-specified knots). ! The last dimension of `fcn` is the output index. ! This is a wrapper for [[db1ink]] (multiple-output version). pure subroutine initialize_1d_specify_knots_multi(me,x,fcn,kx,tx,iflag,extrap) implicit none class(bspline_1d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:,:),intent(in) :: fcn !! `(nx,nf)` array of function values to interpolate !! (the last dimension is the output index) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. integer(ip),intent(out) :: iflag !! status flag (see [[db1ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(1) :: n !! number of points in each direction integer(ip) :: nf !! number of outputs call me%destroy() n(1) = size(x,kind=ip) nf = size(fcn,2,kind=ip) call check_knot_vectors_sizes(nx=n(1),kx=kx,tx=tx,iflag=iflag) if (iflag/=0_ip) then me%iflag = iflag return end if me%nx = n(1) me%kx = kx allocate(me%tx(n(1)+kx)) allocate(me%bcoef(n(1),nf)) allocate(me%work_val_1(3_ip*kx)) me%tx = tx call db1ink(x,n(1),fcn,kx,1_ip,me%tx,me%bcoef,iflag) if (iflag==0_ip) then me%nf = nf call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_1d_specify_knots_multi !***************************************************************************************** !***************************************************************************************** !> ! Evaluate all the outputs of a [[bspline_1d]] with multiple outputs. ! This is a wrapper for [[db1val]] (multiple-output version). pure subroutine evaluate_1d_multi(me,xval,idx,f,iflag) implicit none class(bspline_1d),intent(inout) :: me !! the spline real(wp),intent(in) :: xval !! \(x\) coordinate of evaluation point. integer(ip),intent(in) :: idx !! \(x\) derivative of piecewise polynomial to evaluate. real(wp),dimension(:),intent(out) :: f !! `(nf)` interpolated values of the outputs !! (see [[number_of_outputs]]) integer(ip),intent(out) :: iflag !! status flag (see [[db1val]]) if (me%initialized) then call db1val(xval,idx,me%tx,& me%nx,me%kx,& me%bcoef,f,iflag,me%inbvx,extrap=me%extrap) else f = 0.0_wp iflag = 1_ip end if me%iflag = iflag end subroutine evaluate_1d_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_1d]] type with multiple outputs directly from its ! knots and b-spline coefficients. The last dimension of `bcoef` is the ! output index. pure subroutine initialize_from_coefficients_1d_multi(me,tx,bcoef,kx,iflag,extrap) implicit none class(bspline_1d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. real(wp),dimension(:,:),intent(in) :: bcoef !! `(nx,nf)` array of b-spline coefficients integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(out) :: iflag !! status flag (see [[initialize_from_coefficients_1d]]), !! or 6005 if `bcoef` has no outputs logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) if (size(bcoef,2,kind=ip)<1_ip) then ! no outputs call me%destroy() iflag = 6005_ip me%iflag = iflag return end if call me%initialize_from_coefficients(tx,bcoef(:,1),kx,iflag,extrap) if (iflag==0_ip) then me%bcoef = bcoef me%nf = size(bcoef,2,kind=ip) end if end subroutine initialize_from_coefficients_1d_multi !***************************************************************************************** !***************************************************************************************** !> ! Returns the b-spline coefficients of a [[bspline_1d]] with any number of ! outputs. The last dimension of `bcoef` is the output index. pure subroutine get_coefficients_1d_multi(me,bcoef,iflag) implicit none class(bspline_1d),intent(inout) :: me !! the spline real(wp),dimension(:,:),allocatable,intent(out) :: bcoef !! `(nx,nf)` array of b-spline coefficients integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized) then bcoef = me%bcoef iflag = 0_ip else allocate(bcoef(0,0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_coefficients_1d_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_2d]] type with multiple outputs: `nf` functions of ! the same inputs, on the same grid. ! The last dimension of `fcn` is the output index. ! This is a wrapper for [[db2ink]] (multiple-output version). pure subroutine initialize_2d_auto_knots_multi(me,x,y,fcn,kx,ky,iflag,extrap) implicit none class(bspline_2d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:),intent(in) :: fcn !! `(nx,ny,nf)` array of function values to interpolate !! (the last dimension is the output index) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) integer(ip),intent(out) :: iflag !! status flag (see [[db2ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(2) :: n !! number of points in each direction integer(ip) :: nf !! number of outputs call me%destroy() n(1) = size(x,kind=ip) n(2) = size(y,kind=ip) nf = size(fcn,3,kind=ip) me%nx = n(1) me%ny = n(2) me%kx = kx me%ky = ky allocate(me%tx(n(1)+kx)) allocate(me%ty(n(2)+ky)) allocate(me%bcoef(n(1),n(2),nf)) allocate(me%work_val_1(ky)) allocate(me%work_val_2(3_ip*max(kx,ky))) call db2ink(x,n(1),y,n(2),fcn,kx,ky,0_ip,me%tx,me%ty,me%bcoef,iflag) if (iflag==0_ip) then me%nf = nf call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_2d_auto_knots_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_2d]] type with multiple outputs: `nf` functions of ! the same inputs, on the same grid (with user-specified knots). ! The last dimension of `fcn` is the output index. ! This is a wrapper for [[db2ink]] (multiple-output version). pure subroutine initialize_2d_specify_knots_multi(me,x,y,fcn,kx,ky,tx,ty,iflag,extrap) implicit none class(bspline_2d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:),intent(in) :: fcn !! `(nx,ny,nf)` array of function values to interpolate !! (the last dimension is the output index) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. integer(ip),intent(out) :: iflag !! status flag (see [[db2ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(2) :: n !! number of points in each direction integer(ip) :: nf !! number of outputs call me%destroy() n(1) = size(x,kind=ip) n(2) = size(y,kind=ip) nf = size(fcn,3,kind=ip) call check_knot_vectors_sizes(nx=n(1),kx=kx,tx=tx,ny=n(2),ky=ky,ty=ty,iflag=iflag) if (iflag/=0_ip) then me%iflag = iflag return end if me%nx = n(1) me%ny = n(2) me%kx = kx me%ky = ky allocate(me%tx(n(1)+kx)) allocate(me%ty(n(2)+ky)) allocate(me%bcoef(n(1),n(2),nf)) allocate(me%work_val_1(ky)) allocate(me%work_val_2(3_ip*max(kx,ky))) me%tx = tx me%ty = ty call db2ink(x,n(1),y,n(2),fcn,kx,ky,1_ip,me%tx,me%ty,me%bcoef,iflag) if (iflag==0_ip) then me%nf = nf call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_2d_specify_knots_multi !***************************************************************************************** !***************************************************************************************** !> ! Evaluate all the outputs of a [[bspline_2d]] with multiple outputs. ! This is a wrapper for [[db2val]] (multiple-output version). pure subroutine evaluate_2d_multi(me,xval,yval,idx,idy,f,iflag) implicit none class(bspline_2d),intent(inout) :: me !! the spline real(wp),intent(in) :: xval !! \(x\) coordinate of evaluation point. real(wp),intent(in) :: yval !! \(y\) coordinate of evaluation point. integer(ip),intent(in) :: idx !! \(x\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idy !! \(y\) derivative of piecewise polynomial to evaluate. real(wp),dimension(:),intent(out) :: f !! `(nf)` interpolated values of the outputs !! (see [[number_of_outputs]]) integer(ip),intent(out) :: iflag !! status flag (see [[db2val]]) if (me%initialized) then call db2val(xval,yval,idx,idy,me%tx,me%ty,& me%nx,me%ny,me%kx,me%ky,& me%bcoef,f,iflag,me%inbvx,me%iloy,& me%work_val_1,extrap=me%extrap) else f = 0.0_wp iflag = 1_ip end if me%iflag = iflag end subroutine evaluate_2d_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_2d]] type with multiple outputs directly from its ! knots and b-spline coefficients. The last dimension of `bcoef` is the ! output index. pure subroutine initialize_from_coefficients_2d_multi(me,tx,ty,bcoef,kx,ky,iflag,extrap) implicit none class(bspline_2d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. real(wp),dimension(:,:,:),intent(in) :: bcoef !! `(nx,ny,nf)` array of b-spline coefficients integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) integer(ip),intent(out) :: iflag !! status flag (see [[initialize_from_coefficients_2d]]), !! or 6005 if `bcoef` has no outputs logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) if (size(bcoef,3,kind=ip)<1_ip) then ! no outputs call me%destroy() iflag = 6005_ip me%iflag = iflag return end if call me%initialize_from_coefficients(tx,ty,bcoef(:,:,1),kx,ky,iflag,extrap) if (iflag==0_ip) then me%bcoef = bcoef me%nf = size(bcoef,3,kind=ip) end if end subroutine initialize_from_coefficients_2d_multi !***************************************************************************************** !***************************************************************************************** !> ! Returns the b-spline coefficients of a [[bspline_2d]] with any number of ! outputs. The last dimension of `bcoef` is the output index. pure subroutine get_coefficients_2d_multi(me,bcoef,iflag) implicit none class(bspline_2d),intent(inout) :: me !! the spline real(wp),dimension(:,:,:),allocatable,intent(out) :: bcoef !! `(nx,ny,nf)` array of b-spline coefficients integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized) then bcoef = me%bcoef iflag = 0_ip else allocate(bcoef(0,0,0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_coefficients_2d_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_3d]] type with multiple outputs: `nf` functions of ! the same inputs, on the same grid. ! The last dimension of `fcn` is the output index. ! This is a wrapper for [[db3ink]] (multiple-output version). pure subroutine initialize_3d_auto_knots_multi(me,x,y,z,fcn,kx,ky,kz,iflag,extrap) implicit none class(bspline_3d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nf)` array of function values to interpolate !! (the last dimension is the output index) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) integer(ip),intent(out) :: iflag !! status flag (see [[db3ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(3) :: n !! number of points in each direction integer(ip) :: nf !! number of outputs call me%destroy() n(1) = size(x,kind=ip) n(2) = size(y,kind=ip) n(3) = size(z,kind=ip) nf = size(fcn,4,kind=ip) me%nx = n(1) me%ny = n(2) me%nz = n(3) me%kx = kx me%ky = ky me%kz = kz allocate(me%tx(n(1)+kx)) allocate(me%ty(n(2)+ky)) allocate(me%tz(n(3)+kz)) allocate(me%bcoef(n(1),n(2),n(3),nf)) allocate(me%work_val_1(ky,kz)) allocate(me%work_val_2(kz)) allocate(me%work_val_3(3_ip*max(kx,ky,kz))) call db3ink(x,n(1),y,n(2),z,n(3),fcn,kx,ky,kz,0_ip,me%tx,me%ty,me%tz,me%bcoef,iflag) if (iflag==0_ip) then me%nf = nf call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_3d_auto_knots_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_3d]] type with multiple outputs: `nf` functions of ! the same inputs, on the same grid (with user-specified knots). ! The last dimension of `fcn` is the output index. ! This is a wrapper for [[db3ink]] (multiple-output version). pure subroutine initialize_3d_specify_knots_multi(me,x,y,z,fcn,kx,ky,kz,tx,ty,tz,iflag,extrap) implicit none class(bspline_3d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nf)` array of function values to interpolate !! (the last dimension is the output index) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction. integer(ip),intent(out) :: iflag !! status flag (see [[db3ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(3) :: n !! number of points in each direction integer(ip) :: nf !! number of outputs call me%destroy() n(1) = size(x,kind=ip) n(2) = size(y,kind=ip) n(3) = size(z,kind=ip) nf = size(fcn,4,kind=ip) call check_knot_vectors_sizes(nx=n(1),kx=kx,tx=tx,ny=n(2),ky=ky,ty=ty,nz=n(3),kz=kz,tz=tz,iflag=iflag) if (iflag/=0_ip) then me%iflag = iflag return end if me%nx = n(1) me%ny = n(2) me%nz = n(3) me%kx = kx me%ky = ky me%kz = kz allocate(me%tx(n(1)+kx)) allocate(me%ty(n(2)+ky)) allocate(me%tz(n(3)+kz)) allocate(me%bcoef(n(1),n(2),n(3),nf)) allocate(me%work_val_1(ky,kz)) allocate(me%work_val_2(kz)) allocate(me%work_val_3(3_ip*max(kx,ky,kz))) me%tx = tx me%ty = ty me%tz = tz call db3ink(x,n(1),y,n(2),z,n(3),fcn,kx,ky,kz,1_ip,me%tx,me%ty,me%tz,me%bcoef,iflag) if (iflag==0_ip) then me%nf = nf call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_3d_specify_knots_multi !***************************************************************************************** !***************************************************************************************** !> ! Evaluate all the outputs of a [[bspline_3d]] with multiple outputs. ! This is a wrapper for [[db3val]] (multiple-output version). pure subroutine evaluate_3d_multi(me,xval,yval,zval,idx,idy,idz,f,iflag) implicit none class(bspline_3d),intent(inout) :: me !! the spline real(wp),intent(in) :: xval !! \(x\) coordinate of evaluation point. real(wp),intent(in) :: yval !! \(y\) coordinate of evaluation point. real(wp),intent(in) :: zval !! \(z\) coordinate of evaluation point. integer(ip),intent(in) :: idx !! \(x\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idy !! \(y\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idz !! \(z\) derivative of piecewise polynomial to evaluate. real(wp),dimension(:),intent(out) :: f !! `(nf)` interpolated values of the outputs !! (see [[number_of_outputs]]) integer(ip),intent(out) :: iflag !! status flag (see [[db3val]]) if (me%initialized) then call db3val(xval,yval,zval,idx,idy,idz,me%tx,me%ty,me%tz,& me%nx,me%ny,me%nz,me%kx,me%ky,me%kz,& me%bcoef,f,iflag,me%inbvx,me%iloy,me%iloz,& me%work_val_1,me%work_val_2,extrap=me%extrap) else f = 0.0_wp iflag = 1_ip end if me%iflag = iflag end subroutine evaluate_3d_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_3d]] type with multiple outputs directly from its ! knots and b-spline coefficients. The last dimension of `bcoef` is the ! output index. pure subroutine initialize_from_coefficients_3d_multi(me,tx,ty,tz,bcoef,kx,ky,kz,iflag,extrap) implicit none class(bspline_3d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction. real(wp),dimension(:,:,:,:),intent(in) :: bcoef !! `(nx,ny,nz,nf)` array of b-spline coefficients integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) integer(ip),intent(out) :: iflag !! status flag (see [[initialize_from_coefficients_3d]]), !! or 6005 if `bcoef` has no outputs logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) if (size(bcoef,4,kind=ip)<1_ip) then ! no outputs call me%destroy() iflag = 6005_ip me%iflag = iflag return end if call me%initialize_from_coefficients(tx,ty,tz,bcoef(:,:,:,1),kx,ky,kz,iflag,extrap) if (iflag==0_ip) then me%bcoef = bcoef me%nf = size(bcoef,4,kind=ip) end if end subroutine initialize_from_coefficients_3d_multi !***************************************************************************************** !***************************************************************************************** !> ! Returns the b-spline coefficients of a [[bspline_3d]] with any number of ! outputs. The last dimension of `bcoef` is the output index. pure subroutine get_coefficients_3d_multi(me,bcoef,iflag) implicit none class(bspline_3d),intent(inout) :: me !! the spline real(wp),dimension(:,:,:,:),allocatable,intent(out) :: bcoef !! `(nx,ny,nz,nf)` array of b-spline coefficients integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized) then bcoef = me%bcoef iflag = 0_ip else allocate(bcoef(0,0,0,0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_coefficients_3d_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_4d]] type with multiple outputs: `nf` functions of ! the same inputs, on the same grid. ! The last dimension of `fcn` is the output index. ! This is a wrapper for [[db4ink]] (multiple-output version). pure subroutine initialize_4d_auto_knots_multi(me,x,y,z,q,fcn,kx,ky,kz,kq,iflag,extrap) implicit none class(bspline_4d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nf)` array of function values to interpolate !! (the last dimension is the output index) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) integer(ip),intent(out) :: iflag !! status flag (see [[db4ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(4) :: n !! number of points in each direction integer(ip) :: nf !! number of outputs call me%destroy() n(1) = size(x,kind=ip) n(2) = size(y,kind=ip) n(3) = size(z,kind=ip) n(4) = size(q,kind=ip) nf = size(fcn,5,kind=ip) me%nx = n(1) me%ny = n(2) me%nz = n(3) me%nq = n(4) me%kx = kx me%ky = ky me%kz = kz me%kq = kq allocate(me%tx(n(1)+kx)) allocate(me%ty(n(2)+ky)) allocate(me%tz(n(3)+kz)) allocate(me%tq(n(4)+kq)) allocate(me%bcoef(n(1),n(2),n(3),n(4),nf)) allocate(me%work_val_1(ky,kz,kq)) allocate(me%work_val_2(kz,kq)) allocate(me%work_val_3(kq)) allocate(me%work_val_4(3_ip*max(kx,ky,kz,kq))) call db4ink(x,n(1),y,n(2),z,n(3),q,n(4),fcn,kx,ky,kz,kq,0_ip,me%tx,me%ty,me%tz,me%tq,me%bcoef,iflag) if (iflag==0_ip) then me%nf = nf call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_4d_auto_knots_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_4d]] type with multiple outputs: `nf` functions of ! the same inputs, on the same grid (with user-specified knots). ! The last dimension of `fcn` is the output index. ! This is a wrapper for [[db4ink]] (multiple-output version). pure subroutine initialize_4d_specify_knots_multi(me,x,y,z,q,fcn,kx,ky,kz,kq,tx,ty,tz,tq,iflag,extrap) implicit none class(bspline_4d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nf)` array of function values to interpolate !! (the last dimension is the output index) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction. integer(ip),intent(out) :: iflag !! status flag (see [[db4ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(4) :: n !! number of points in each direction integer(ip) :: nf !! number of outputs call me%destroy() n(1) = size(x,kind=ip) n(2) = size(y,kind=ip) n(3) = size(z,kind=ip) n(4) = size(q,kind=ip) nf = size(fcn,5,kind=ip) call check_knot_vectors_sizes(nx=n(1),kx=kx,tx=tx,ny=n(2),ky=ky,ty=ty,nz=n(3),kz=kz,tz=tz,nq=n(4),kq=kq,tq=tq,iflag=iflag) if (iflag/=0_ip) then me%iflag = iflag return end if me%nx = n(1) me%ny = n(2) me%nz = n(3) me%nq = n(4) me%kx = kx me%ky = ky me%kz = kz me%kq = kq allocate(me%tx(n(1)+kx)) allocate(me%ty(n(2)+ky)) allocate(me%tz(n(3)+kz)) allocate(me%tq(n(4)+kq)) allocate(me%bcoef(n(1),n(2),n(3),n(4),nf)) allocate(me%work_val_1(ky,kz,kq)) allocate(me%work_val_2(kz,kq)) allocate(me%work_val_3(kq)) allocate(me%work_val_4(3_ip*max(kx,ky,kz,kq))) me%tx = tx me%ty = ty me%tz = tz me%tq = tq call db4ink(x,n(1),y,n(2),z,n(3),q,n(4),fcn,kx,ky,kz,kq,1_ip,me%tx,me%ty,me%tz,me%tq,me%bcoef,iflag) if (iflag==0_ip) then me%nf = nf call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_4d_specify_knots_multi !***************************************************************************************** !***************************************************************************************** !> ! Evaluate all the outputs of a [[bspline_4d]] with multiple outputs. ! This is a wrapper for [[db4val]] (multiple-output version). pure subroutine evaluate_4d_multi(me,xval,yval,zval,qval,idx,idy,idz,idq,f,iflag) implicit none class(bspline_4d),intent(inout) :: me !! the spline real(wp),intent(in) :: xval !! \(x\) coordinate of evaluation point. real(wp),intent(in) :: yval !! \(y\) coordinate of evaluation point. real(wp),intent(in) :: zval !! \(z\) coordinate of evaluation point. real(wp),intent(in) :: qval !! \(q\) coordinate of evaluation point. integer(ip),intent(in) :: idx !! \(x\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idy !! \(y\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idz !! \(z\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idq !! \(q\) derivative of piecewise polynomial to evaluate. real(wp),dimension(:),intent(out) :: f !! `(nf)` interpolated values of the outputs !! (see [[number_of_outputs]]) integer(ip),intent(out) :: iflag !! status flag (see [[db4val]]) if (me%initialized) then call db4val(xval,yval,zval,qval,idx,idy,idz,idq,me%tx,me%ty,me%tz,me%tq,& me%nx,me%ny,me%nz,me%nq,me%kx,me%ky,me%kz,me%kq,& me%bcoef,f,iflag,me%inbvx,me%iloy,me%iloz,me%iloq,& me%work_val_1,me%work_val_2,me%work_val_3,extrap=me%extrap) else f = 0.0_wp iflag = 1_ip end if me%iflag = iflag end subroutine evaluate_4d_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_4d]] type with multiple outputs directly from its ! knots and b-spline coefficients. The last dimension of `bcoef` is the ! output index. pure subroutine initialize_from_coefficients_4d_multi(me,tx,ty,tz,tq,bcoef,kx,ky,kz,kq,iflag,extrap) implicit none class(bspline_4d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction. real(wp),dimension(:,:,:,:,:),intent(in) :: bcoef !! `(nx,ny,nz,nq,nf)` array of b-spline coefficients integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) integer(ip),intent(out) :: iflag !! status flag (see [[initialize_from_coefficients_4d]]), !! or 6005 if `bcoef` has no outputs logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) if (size(bcoef,5,kind=ip)<1_ip) then ! no outputs call me%destroy() iflag = 6005_ip me%iflag = iflag return end if call me%initialize_from_coefficients(tx,ty,tz,tq,bcoef(:,:,:,:,1),kx,ky,kz,kq,iflag,extrap) if (iflag==0_ip) then me%bcoef = bcoef me%nf = size(bcoef,5,kind=ip) end if end subroutine initialize_from_coefficients_4d_multi !***************************************************************************************** !***************************************************************************************** !> ! Returns the b-spline coefficients of a [[bspline_4d]] with any number of ! outputs. The last dimension of `bcoef` is the output index. pure subroutine get_coefficients_4d_multi(me,bcoef,iflag) implicit none class(bspline_4d),intent(inout) :: me !! the spline real(wp),dimension(:,:,:,:,:),allocatable,intent(out) :: bcoef !! `(nx,ny,nz,nq,nf)` array of b-spline coefficients integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized) then bcoef = me%bcoef iflag = 0_ip else allocate(bcoef(0,0,0,0,0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_coefficients_4d_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_5d]] type with multiple outputs: `nf` functions of ! the same inputs, on the same grid. ! The last dimension of `fcn` is the output index. ! This is a wrapper for [[db5ink]] (multiple-output version). pure subroutine initialize_5d_auto_knots_multi(me,x,y,z,q,r,fcn,kx,ky,kz,kq,kr,iflag,extrap) implicit none class(bspline_5d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: r !! `(nr)` array of \(r\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr,nf)` array of function values to interpolate !! (the last dimension is the output index) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) integer(ip),intent(out) :: iflag !! status flag (see [[db5ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(5) :: n !! number of points in each direction integer(ip) :: nf !! number of outputs call me%destroy() n(1) = size(x,kind=ip) n(2) = size(y,kind=ip) n(3) = size(z,kind=ip) n(4) = size(q,kind=ip) n(5) = size(r,kind=ip) nf = size(fcn,6,kind=ip) me%nx = n(1) me%ny = n(2) me%nz = n(3) me%nq = n(4) me%nr = n(5) me%kx = kx me%ky = ky me%kz = kz me%kq = kq me%kr = kr allocate(me%tx(n(1)+kx)) allocate(me%ty(n(2)+ky)) allocate(me%tz(n(3)+kz)) allocate(me%tq(n(4)+kq)) allocate(me%tr(n(5)+kr)) allocate(me%bcoef(n(1),n(2),n(3),n(4),n(5),nf)) allocate(me%work_val_1(ky,kz,kq,kr)) allocate(me%work_val_2(kz,kq,kr)) allocate(me%work_val_3(kq,kr)) allocate(me%work_val_4(kr)) allocate(me%work_val_5(3_ip*max(kx,ky,kz,kq,kr))) call db5ink(x,n(1),y,n(2),z,n(3),q,n(4),r,n(5),fcn,kx,ky,kz,kq,kr,0_ip,me%tx,me%ty,me%tz,me%tq,me%tr,me%bcoef,iflag) if (iflag==0_ip) then me%nf = nf call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_5d_auto_knots_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_5d]] type with multiple outputs: `nf` functions of ! the same inputs, on the same grid (with user-specified knots). ! The last dimension of `fcn` is the output index. ! This is a wrapper for [[db5ink]] (multiple-output version). pure subroutine initialize_5d_specify_knots_multi(me,x,y,z,q,r,fcn,kx,ky,kz,kq,kr,tx,ty,tz,tq,tr,iflag,extrap) implicit none class(bspline_5d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: r !! `(nr)` array of \(r\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr,nf)` array of function values to interpolate !! (the last dimension is the output index) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction. real(wp),dimension(:),intent(in) :: tr !! The `(nr+kr)` knots in the \(r\) direction. integer(ip),intent(out) :: iflag !! status flag (see [[db5ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(5) :: n !! number of points in each direction integer(ip) :: nf !! number of outputs call me%destroy() n(1) = size(x,kind=ip) n(2) = size(y,kind=ip) n(3) = size(z,kind=ip) n(4) = size(q,kind=ip) n(5) = size(r,kind=ip) nf = size(fcn,6,kind=ip) call check_knot_vectors_sizes(nx=n(1),kx=kx,tx=tx,ny=n(2),ky=ky,ty=ty,nz=n(3),kz=kz,tz=tz,nq=n(4),kq=kq,tq=tq,nr=n(5),kr=kr,tr=tr,iflag=iflag) if (iflag/=0_ip) then me%iflag = iflag return end if me%nx = n(1) me%ny = n(2) me%nz = n(3) me%nq = n(4) me%nr = n(5) me%kx = kx me%ky = ky me%kz = kz me%kq = kq me%kr = kr allocate(me%tx(n(1)+kx)) allocate(me%ty(n(2)+ky)) allocate(me%tz(n(3)+kz)) allocate(me%tq(n(4)+kq)) allocate(me%tr(n(5)+kr)) allocate(me%bcoef(n(1),n(2),n(3),n(4),n(5),nf)) allocate(me%work_val_1(ky,kz,kq,kr)) allocate(me%work_val_2(kz,kq,kr)) allocate(me%work_val_3(kq,kr)) allocate(me%work_val_4(kr)) allocate(me%work_val_5(3_ip*max(kx,ky,kz,kq,kr))) me%tx = tx me%ty = ty me%tz = tz me%tq = tq me%tr = tr call db5ink(x,n(1),y,n(2),z,n(3),q,n(4),r,n(5),fcn,kx,ky,kz,kq,kr,1_ip,me%tx,me%ty,me%tz,me%tq,me%tr,me%bcoef,iflag) if (iflag==0_ip) then me%nf = nf call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_5d_specify_knots_multi !***************************************************************************************** !***************************************************************************************** !> ! Evaluate all the outputs of a [[bspline_5d]] with multiple outputs. ! This is a wrapper for [[db5val]] (multiple-output version). pure subroutine evaluate_5d_multi(me,xval,yval,zval,qval,rval,idx,idy,idz,idq,idr,f,iflag) implicit none class(bspline_5d),intent(inout) :: me !! the spline real(wp),intent(in) :: xval !! \(x\) coordinate of evaluation point. real(wp),intent(in) :: yval !! \(y\) coordinate of evaluation point. real(wp),intent(in) :: zval !! \(z\) coordinate of evaluation point. real(wp),intent(in) :: qval !! \(q\) coordinate of evaluation point. real(wp),intent(in) :: rval !! \(r\) coordinate of evaluation point. integer(ip),intent(in) :: idx !! \(x\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idy !! \(y\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idz !! \(z\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idq !! \(q\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idr !! \(r\) derivative of piecewise polynomial to evaluate. real(wp),dimension(:),intent(out) :: f !! `(nf)` interpolated values of the outputs !! (see [[number_of_outputs]]) integer(ip),intent(out) :: iflag !! status flag (see [[db5val]]) if (me%initialized) then call db5val(xval,yval,zval,qval,rval,idx,idy,idz,idq,idr,me%tx,me%ty,me%tz,me%tq,me%tr,& me%nx,me%ny,me%nz,me%nq,me%nr,me%kx,me%ky,me%kz,me%kq,me%kr,& me%bcoef,f,iflag,me%inbvx,me%iloy,me%iloz,me%iloq,me%ilor,& me%work_val_1,me%work_val_2,me%work_val_3,me%work_val_4,extrap=me%extrap) else f = 0.0_wp iflag = 1_ip end if me%iflag = iflag end subroutine evaluate_5d_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_5d]] type with multiple outputs directly from its ! knots and b-spline coefficients. The last dimension of `bcoef` is the ! output index. pure subroutine initialize_from_coefficients_5d_multi(me,tx,ty,tz,tq,tr,bcoef,kx,ky,kz,kq,kr,iflag,extrap) implicit none class(bspline_5d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction. real(wp),dimension(:),intent(in) :: tr !! The `(nr+kr)` knots in the \(r\) direction. real(wp),dimension(:,:,:,:,:,:),intent(in) :: bcoef !! `(nx,ny,nz,nq,nr,nf)` array of b-spline coefficients integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) integer(ip),intent(out) :: iflag !! status flag (see [[initialize_from_coefficients_5d]]), !! or 6005 if `bcoef` has no outputs logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) if (size(bcoef,6,kind=ip)<1_ip) then ! no outputs call me%destroy() iflag = 6005_ip me%iflag = iflag return end if call me%initialize_from_coefficients(tx,ty,tz,tq,tr,bcoef(:,:,:,:,:,1),kx,ky,kz,kq,kr,iflag,extrap) if (iflag==0_ip) then me%bcoef = bcoef me%nf = size(bcoef,6,kind=ip) end if end subroutine initialize_from_coefficients_5d_multi !***************************************************************************************** !***************************************************************************************** !> ! Returns the b-spline coefficients of a [[bspline_5d]] with any number of ! outputs. The last dimension of `bcoef` is the output index. pure subroutine get_coefficients_5d_multi(me,bcoef,iflag) implicit none class(bspline_5d),intent(inout) :: me !! the spline real(wp),dimension(:,:,:,:,:,:),allocatable,intent(out) :: bcoef !! `(nx,ny,nz,nq,nr,nf)` array of b-spline coefficients integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized) then bcoef = me%bcoef iflag = 0_ip else allocate(bcoef(0,0,0,0,0,0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_coefficients_5d_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_6d]] type with multiple outputs: `nf` functions of ! the same inputs, on the same grid. ! The last dimension of `fcn` is the output index. ! This is a wrapper for [[db6ink]] (multiple-output version). pure subroutine initialize_6d_auto_knots_multi(me,x,y,z,q,r,s,fcn,kx,ky,kz,kq,kr,ks,iflag,extrap) implicit none class(bspline_6d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: r !! `(nr)` array of \(r\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: s !! `(ns)` array of \(s\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr,ns,nf)` array of function values to interpolate !! (the last dimension is the output index) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) integer(ip),intent(in) :: ks !! The order of spline pieces in \(s\) integer(ip),intent(out) :: iflag !! status flag (see [[db6ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(6) :: n !! number of points in each direction integer(ip) :: nf !! number of outputs call me%destroy() n(1) = size(x,kind=ip) n(2) = size(y,kind=ip) n(3) = size(z,kind=ip) n(4) = size(q,kind=ip) n(5) = size(r,kind=ip) n(6) = size(s,kind=ip) nf = size(fcn,7,kind=ip) me%nx = n(1) me%ny = n(2) me%nz = n(3) me%nq = n(4) me%nr = n(5) me%ns = n(6) me%kx = kx me%ky = ky me%kz = kz me%kq = kq me%kr = kr me%ks = ks allocate(me%tx(n(1)+kx)) allocate(me%ty(n(2)+ky)) allocate(me%tz(n(3)+kz)) allocate(me%tq(n(4)+kq)) allocate(me%tr(n(5)+kr)) allocate(me%ts(n(6)+ks)) allocate(me%bcoef(n(1),n(2),n(3),n(4),n(5),n(6),nf)) allocate(me%work_val_1(ky,kz,kq,kr,ks)) allocate(me%work_val_2(kz,kq,kr,ks)) allocate(me%work_val_3(kq,kr,ks)) allocate(me%work_val_4(kr,ks)) allocate(me%work_val_5(ks)) allocate(me%work_val_6(3_ip*max(kx,ky,kz,kq,kr,ks))) call db6ink(x,n(1),y,n(2),z,n(3),q,n(4),r,n(5),s,n(6),fcn,kx,ky,kz,kq,kr,ks,0_ip,me%tx,me%ty,me%tz,me%tq,me%tr,me%ts,me%bcoef,iflag) if (iflag==0_ip) then me%nf = nf call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_6d_auto_knots_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_6d]] type with multiple outputs: `nf` functions of ! the same inputs, on the same grid (with user-specified knots). ! The last dimension of `fcn` is the output index. ! This is a wrapper for [[db6ink]] (multiple-output version). pure subroutine initialize_6d_specify_knots_multi(me,x,y,z,q,r,s,fcn,kx,ky,kz,kq,kr,ks,tx,ty,tz,tq,tr,ts,iflag,extrap) implicit none class(bspline_6d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: x !! `(nx)` array of \(x\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: y !! `(ny)` array of \(y\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: z !! `(nz)` array of \(z\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: q !! `(nq)` array of \(q\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: r !! `(nr)` array of \(r\) abcissae. Must be strictly increasing. real(wp),dimension(:),intent(in) :: s !! `(ns)` array of \(s\) abcissae. Must be strictly increasing. real(wp),dimension(:,:,:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr,ns,nf)` array of function values to interpolate !! (the last dimension is the output index) integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) integer(ip),intent(in) :: ks !! The order of spline pieces in \(s\) real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction. real(wp),dimension(:),intent(in) :: tr !! The `(nr+kr)` knots in the \(r\) direction. real(wp),dimension(:),intent(in) :: ts !! The `(ns+ks)` knots in the \(s\) direction. integer(ip),intent(out) :: iflag !! status flag (see [[db6ink]]) logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) integer(ip),dimension(6) :: n !! number of points in each direction integer(ip) :: nf !! number of outputs call me%destroy() n(1) = size(x,kind=ip) n(2) = size(y,kind=ip) n(3) = size(z,kind=ip) n(4) = size(q,kind=ip) n(5) = size(r,kind=ip) n(6) = size(s,kind=ip) nf = size(fcn,7,kind=ip) call check_knot_vectors_sizes(nx=n(1),kx=kx,tx=tx,ny=n(2),ky=ky,ty=ty,nz=n(3),kz=kz,tz=tz,nq=n(4),kq=kq,tq=tq,nr=n(5),kr=kr,tr=tr,ns=n(6),ks=ks,ts=ts,iflag=iflag) if (iflag/=0_ip) then me%iflag = iflag return end if me%nx = n(1) me%ny = n(2) me%nz = n(3) me%nq = n(4) me%nr = n(5) me%ns = n(6) me%kx = kx me%ky = ky me%kz = kz me%kq = kq me%kr = kr me%ks = ks allocate(me%tx(n(1)+kx)) allocate(me%ty(n(2)+ky)) allocate(me%tz(n(3)+kz)) allocate(me%tq(n(4)+kq)) allocate(me%tr(n(5)+kr)) allocate(me%ts(n(6)+ks)) allocate(me%bcoef(n(1),n(2),n(3),n(4),n(5),n(6),nf)) allocate(me%work_val_1(ky,kz,kq,kr,ks)) allocate(me%work_val_2(kz,kq,kr,ks)) allocate(me%work_val_3(kq,kr,ks)) allocate(me%work_val_4(kr,ks)) allocate(me%work_val_5(ks)) allocate(me%work_val_6(3_ip*max(kx,ky,kz,kq,kr,ks))) me%tx = tx me%ty = ty me%tz = tz me%tq = tq me%tr = tr me%ts = ts call db6ink(x,n(1),y,n(2),z,n(3),q,n(4),r,n(5),s,n(6),fcn,kx,ky,kz,kq,kr,ks,1_ip,me%tx,me%ty,me%tz,me%tq,me%tr,me%ts,me%bcoef,iflag) if (iflag==0_ip) then me%nf = nf call me%set_extrap_flag(extrap) end if me%initialized = iflag==0_ip me%iflag = iflag end subroutine initialize_6d_specify_knots_multi !***************************************************************************************** !***************************************************************************************** !> ! Evaluate all the outputs of a [[bspline_6d]] with multiple outputs. ! This is a wrapper for [[db6val]] (multiple-output version). pure subroutine evaluate_6d_multi(me,xval,yval,zval,qval,rval,sval,idx,idy,idz,idq,idr,ids,f,iflag) implicit none class(bspline_6d),intent(inout) :: me !! the spline real(wp),intent(in) :: xval !! \(x\) coordinate of evaluation point. real(wp),intent(in) :: yval !! \(y\) coordinate of evaluation point. real(wp),intent(in) :: zval !! \(z\) coordinate of evaluation point. real(wp),intent(in) :: qval !! \(q\) coordinate of evaluation point. real(wp),intent(in) :: rval !! \(r\) coordinate of evaluation point. real(wp),intent(in) :: sval !! \(s\) coordinate of evaluation point. integer(ip),intent(in) :: idx !! \(x\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idy !! \(y\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idz !! \(z\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idq !! \(q\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: idr !! \(r\) derivative of piecewise polynomial to evaluate. integer(ip),intent(in) :: ids !! \(s\) derivative of piecewise polynomial to evaluate. real(wp),dimension(:),intent(out) :: f !! `(nf)` interpolated values of the outputs !! (see [[number_of_outputs]]) integer(ip),intent(out) :: iflag !! status flag (see [[db6val]]) if (me%initialized) then call db6val(xval,yval,zval,qval,rval,sval,idx,idy,idz,idq,idr,ids,me%tx,me%ty,me%tz,me%tq,me%tr,me%ts,& me%nx,me%ny,me%nz,me%nq,me%nr,me%ns,me%kx,me%ky,me%kz,me%kq,me%kr,me%ks,& me%bcoef,f,iflag,me%inbvx,me%iloy,me%iloz,me%iloq,me%ilor,me%ilos,& me%work_val_1,me%work_val_2,me%work_val_3,me%work_val_4,me%work_val_5,extrap=me%extrap) else f = 0.0_wp iflag = 1_ip end if me%iflag = iflag end subroutine evaluate_6d_multi !***************************************************************************************** !***************************************************************************************** !> ! Initialize a [[bspline_6d]] type with multiple outputs directly from its ! knots and b-spline coefficients. The last dimension of `bcoef` is the ! output index. pure subroutine initialize_from_coefficients_6d_multi(me,tx,ty,tz,tq,tr,ts,bcoef,kx,ky,kz,kq,kr,ks,iflag,extrap) implicit none class(bspline_6d),intent(inout) :: me !! the spline real(wp),dimension(:),intent(in) :: tx !! The `(nx+kx)` knots in the \(x\) direction. real(wp),dimension(:),intent(in) :: ty !! The `(ny+ky)` knots in the \(y\) direction. real(wp),dimension(:),intent(in) :: tz !! The `(nz+kz)` knots in the \(z\) direction. real(wp),dimension(:),intent(in) :: tq !! The `(nq+kq)` knots in the \(q\) direction. real(wp),dimension(:),intent(in) :: tr !! The `(nr+kr)` knots in the \(r\) direction. real(wp),dimension(:),intent(in) :: ts !! The `(ns+ks)` knots in the \(s\) direction. real(wp),dimension(:,:,:,:,:,:,:),intent(in) :: bcoef !! `(nx,ny,nz,nq,nr,ns,nf)` array of b-spline coefficients integer(ip),intent(in) :: kx !! The order of spline pieces in \(x\) integer(ip),intent(in) :: ky !! The order of spline pieces in \(y\) integer(ip),intent(in) :: kz !! The order of spline pieces in \(z\) integer(ip),intent(in) :: kq !! The order of spline pieces in \(q\) integer(ip),intent(in) :: kr !! The order of spline pieces in \(r\) integer(ip),intent(in) :: ks !! The order of spline pieces in \(s\) integer(ip),intent(out) :: iflag !! status flag (see [[initialize_from_coefficients_6d]]), !! or 6005 if `bcoef` has no outputs logical,intent(in),optional :: extrap !! if true, then extrapolation is allowed !! (default is false) if (size(bcoef,7,kind=ip)<1_ip) then ! no outputs call me%destroy() iflag = 6005_ip me%iflag = iflag return end if call me%initialize_from_coefficients(tx,ty,tz,tq,tr,ts,bcoef(:,:,:,:,:,:,1),kx,ky,kz,kq,kr,ks,iflag,extrap) if (iflag==0_ip) then me%bcoef = bcoef me%nf = size(bcoef,7,kind=ip) end if end subroutine initialize_from_coefficients_6d_multi !***************************************************************************************** !***************************************************************************************** !> ! Returns the b-spline coefficients of a [[bspline_6d]] with any number of ! outputs. The last dimension of `bcoef` is the output index. pure subroutine get_coefficients_6d_multi(me,bcoef,iflag) implicit none class(bspline_6d),intent(inout) :: me !! the spline real(wp),dimension(:,:,:,:,:,:,:),allocatable,intent(out) :: bcoef !! `(nx,ny,nz,nq,nr,ns,nf)` array of b-spline coefficients integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized) then bcoef = me%bcoef iflag = 0_ip else allocate(bcoef(0,0,0,0,0,0,0)) iflag = 1_ip end if me%iflag = iflag end subroutine get_coefficients_6d_multi !***************************************************************************************** !***************************************************************************************** !> ! The control points of a [[bspline_1d]]: the Greville abscissae ! (see [[greville_abscissae]]) and the b-spline coefficients. ! ! The spline lies within the convex hull of its control points, and ! the control polygon (the control points joined by straight lines) ! approximates the spline. pure subroutine control_points_1d(me,x,y,iflag) implicit none class(bspline_1d),intent(inout) :: me real(wp),dimension(:),allocatable,intent(out) :: x !! `(nx)` locations of the control points !! (the Greville abscissae) real(wp),dimension(:),allocatable,intent(out) :: y !! `(nx)` values of the control points !! (the b-spline coefficients) integer(ip),intent(out) :: iflag !! status flag: 0 if no errors, !! 1 if the class is not initialized if (me%initialized .and. me%nf/=1_ip) then iflag = 6003_ip ! use the multiple-output version allocate(x(0), y(0)) else if (me%initialized) then allocate(x(me%nx)) call greville_abscissae(me%tx,me%nx,me%kx,x) y = me%bcoef(:,1) iflag = 0_ip else allocate(x(0), y(0)) iflag = 1_ip end if me%iflag = iflag end subroutine control_points_1d !***************************************************************************************** !***************************************************************************************** !> ! Check the knots and coefficients for one direction in the ! `initialize_from_coefficients` methods. pure subroutine check_coefficients(n,k,t,iflag) implicit none integer(ip),intent(in) :: n !! number of b-spline coefficients integer(ip),intent(in) :: k !! order of the b-spline real(wp),dimension(:),intent(in) :: t !! knot vector integer(ip),intent(out) :: iflag !! status flag (0 if everything is OK) integer(ip) :: nt !! number of knots nt = size(t,kind=ip) if (k<1_ip) then iflag = 5001_ip else if (n<k) then iflag = 5002_ip else if (nt/=n+k) then iflag = 5003_ip else if (any(t(2_ip:nt)<t(1_ip:nt-1_ip))) then iflag = 5004_ip else if (t(k)>=t(n+1_ip)) then iflag = 5005_ip else iflag = 0_ip end if end subroutine check_coefficients !***************************************************************************************** !***************************************************************************************** end module bspline_oo_module !*****************************************************************************************