db5val_multi Subroutine

private pure subroutine db5val_multi(xval, yval, zval, qval, rval, idx, idy, idz, idq, idr, tx, ty, tz, tq, tr, nx, ny, nz, nq, nr, kx, ky, kz, kq, kr, bcoef, f, iflag, inbvx, iloy, iloz, iloq, ilor, w4, w3, w2, w1, extrap)

Multiple-output version of db5val: evaluates all nf outputs of a spline from db5ink (multiple-output version) at one point. The interval search and the basis functions in each direction are only computed once, for all of the outputs.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: xval

coordinate of evaluation point.
real(kind=wp), intent(in) :: yval

coordinate of evaluation point.
real(kind=wp), intent(in) :: zval

coordinate of evaluation point.
real(kind=wp), intent(in) :: qval

coordinate of evaluation point.
real(kind=wp), intent(in) :: rval

coordinate of evaluation point.
integer(kind=ip), intent(in) :: idx

derivative of piecewise polynomial to evaluate.
integer(kind=ip), intent(in) :: idy

derivative of piecewise polynomial to evaluate.
integer(kind=ip), intent(in) :: idz

derivative of piecewise polynomial to evaluate.
integer(kind=ip), intent(in) :: idq

derivative of piecewise polynomial to evaluate.
integer(kind=ip), intent(in) :: idr

derivative of piecewise polynomial to evaluate.
real(kind=wp), intent(in), dimension(nx+kx) :: tx

knots in the direction.
real(kind=wp), intent(in), dimension(ny+ky) :: ty

knots in the direction.
real(kind=wp), intent(in), dimension(nz+kz) :: tz

knots in the direction.
real(kind=wp), intent(in), dimension(nq+kq) :: tq

knots in the direction.
real(kind=wp), intent(in), dimension(nr+kr) :: tr

knots in the direction.
integer(kind=ip), intent(in) :: nx

the number of interpolation points in .
integer(kind=ip), intent(in) :: ny

the number of interpolation points in .
integer(kind=ip), intent(in) :: nz

the number of interpolation points in .
integer(kind=ip), intent(in) :: nq

the number of interpolation points in .
integer(kind=ip), intent(in) :: nr

the number of interpolation points in .
integer(kind=ip), intent(in) :: kx

order of polynomial pieces in .
integer(kind=ip), intent(in) :: ky

order of polynomial pieces in .
integer(kind=ip), intent(in) :: kz

order of polynomial pieces in .
integer(kind=ip), intent(in) :: kq

order of polynomial pieces in .
integer(kind=ip), intent(in) :: kr

order of polynomial pieces in .
real(kind=wp), intent(in), dimension(:,:,:,:,:,:) :: bcoef

(nx,ny,nz,nq,nr,nf) b-spline coefficients from db5ink (multiple-output version).
real(kind=wp), intent(out), dimension(:) :: f

(nf) interpolated values of the outputs
integer(kind=ip), intent(out) :: iflag

status flag:

  • : no errors
  • 6002: size(f) is not the number of outputs
  • 6004: the leading dimensions of bcoef are not (nx,ny,nz,nq,nr)
  • other : error (see db5val)
integer(kind=ip), intent(inout) :: inbvx

initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.
integer(kind=ip), intent(inout) :: iloy

initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.
integer(kind=ip), intent(inout) :: iloz

initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.
integer(kind=ip), intent(inout) :: iloq

initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.
integer(kind=ip), intent(inout) :: ilor

initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.
real(kind=wp), intent(inout), dimension(ky,kz,kq,kr) :: w4

work array of dimension (ky,kz,kq,kr)
real(kind=wp), intent(inout), dimension(kz,kq,kr) :: w3

work array of dimension (kz,kq,kr)
real(kind=wp), intent(inout), dimension(kq,kr) :: w2

work array of dimension (kq,kr)
real(kind=wp), intent(inout), dimension(kr) :: w1

work array of dimension (kr)
logical, intent(in), optional :: extrap

if extrapolation is allowed (if not present, default is False)

Calls

proc~~db5val_multi~~CallsGraph proc~db5val_multi db5val_multi proc~check_value check_value proc~db5val_multi->proc~check_value proc~dbvalu_weights dbvalu_weights proc~db5val_multi->proc~dbvalu_weights proc~dintrv dintrv proc~dbvalu_weights->proc~dintrv proc~get_temp_x_for_extrap get_temp_x_for_extrap proc~dintrv->proc~get_temp_x_for_extrap

Called by

proc~~db5val_multi~~CalledByGraph proc~db5val_multi db5val_multi interface~db5val db5val interface~db5val->proc~db5val_multi proc~evaluate_5d bspline_5d%evaluate_5d proc~evaluate_5d->interface~db5val proc~evaluate_5d_multi bspline_5d%evaluate_5d_multi proc~evaluate_5d_multi->interface~db5val

Source Code

    pure subroutine db5val_multi(xval,yval,zval,qval,rval,idx,idy,idz,idq,idr,tx,ty,tz,tq,tr,nx,ny,nz,nq,nr,kx,ky,kz,kq,kr,bcoef,f,iflag,inbvx,iloy,iloz,iloq,ilor,w4,w3,w2,w1,extrap)

    implicit none

    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.
    integer(ip),intent(in)              :: nx     !! the number of interpolation points in \(x\).
    integer(ip),intent(in)              :: ny     !! the number of interpolation points in \(y\).
    integer(ip),intent(in)              :: nz     !! the number of interpolation points in \(z\).
    integer(ip),intent(in)              :: nq     !! the number of interpolation points in \(q\).
    integer(ip),intent(in)              :: nr     !! the number of interpolation points in \(r\).
    integer(ip),intent(in)              :: kx     !! order of polynomial pieces in \(x\).
    integer(ip),intent(in)              :: ky     !! order of polynomial pieces in \(y\).
    integer(ip),intent(in)              :: kz     !! order of polynomial pieces in \(z\).
    integer(ip),intent(in)              :: kq     !! order of polynomial pieces in \(q\).
    integer(ip),intent(in)              :: kr     !! order of polynomial pieces in \(r\).
    real(wp),dimension(nx+kx),intent(in) :: tx !! knots in the \(x\) direction.
    real(wp),dimension(ny+ky),intent(in) :: ty !! knots in the \(y\) direction.
    real(wp),dimension(nz+kz),intent(in) :: tz !! knots in the \(z\) direction.
    real(wp),dimension(nq+kq),intent(in) :: tq !! knots in the \(q\) direction.
    real(wp),dimension(nr+kr),intent(in) :: tr !! knots in the \(r\) direction.
    real(wp),dimension(:,:,:,:,:,:),intent(in) :: bcoef  !! `(nx,ny,nz,nq,nr,nf)` b-spline coefficients
                                                  !! from [[db5ink]] (multiple-output version).
    real(wp),dimension(:),intent(out)   :: f      !! `(nf)` interpolated values of the outputs
    integer(ip),intent(out)             :: iflag  !! status flag:
                                                  !!
                                                  !! * \( = 0 \)   : no errors
                                                  !! * 6002: `size(f)` is not the number of outputs
                                                  !! * 6004: the leading dimensions of `bcoef` are not `(nx,ny,nz,nq,nr)`
                                                  !! * other \( \ne 0 \) : error (see [[db5val]])
    integer(ip),intent(inout)           :: inbvx  !! initialization parameter which must be
                                                  !! set to 1 the first time this routine is called,
                                                  !! and must not be changed by the user.
    integer(ip),intent(inout)           :: iloy   !! initialization parameter which must be
                                                  !! set to 1 the first time this routine is called,
                                                  !! and must not be changed by the user.
    integer(ip),intent(inout)           :: iloz   !! initialization parameter which must be
                                                  !! set to 1 the first time this routine is called,
                                                  !! and must not be changed by the user.
    integer(ip),intent(inout)           :: iloq   !! initialization parameter which must be
                                                  !! set to 1 the first time this routine is called,
                                                  !! and must not be changed by the user.
    integer(ip),intent(inout)           :: ilor   !! initialization parameter which must be
                                                  !! set to 1 the first time this routine is called,
                                                  !! and must not be changed by the user.
    logical,intent(in),optional         :: extrap !! if extrapolation is allowed
                                                  !! (if not present, default is False)

    integer(ip) :: leftx  !! interval index in \(x\)
    integer(ip) :: lefty  !! interval index in \(y\)
    integer(ip) :: leftz  !! interval index in \(z\)
    integer(ip) :: leftq  !! interval index in \(q\)
    integer(ip) :: leftr  !! interval index in \(r\)
    integer(ip) :: ix     !! coefficient offset in \(x\)
    integer(ip) :: iy     !! coefficient offset in \(y\)
    integer(ip) :: iz     !! coefficient offset in \(z\)
    integer(ip) :: iq     !! coefficient offset in \(q\)
    integer(ip) :: ir     !! coefficient offset in \(r\)
    real(wp),dimension(kx) :: bx  !! basis function weights in \(x\)
    real(wp),dimension(ky) :: by  !! basis function weights in \(y\)
    real(wp),dimension(kz) :: bz  !! basis function weights in \(z\)
    real(wp),dimension(kq) :: bq  !! basis function weights in \(q\)
    real(wp),dimension(kr) :: br  !! basis function weights in \(r\)
    real(wp),dimension(ky,kz,kq,kr),intent(inout) :: w4  !! work array of dimension `(ky,kz,kq,kr)`
    real(wp),dimension(kz,kq,kr),intent(inout) :: w3  !! work array of dimension `(kz,kq,kr)`
    real(wp),dimension(kq,kr),intent(inout) :: w2  !! work array of dimension `(kq,kr)`
    real(wp),dimension(kr),intent(inout) :: w1  !! work array of dimension `(kr)`
    integer(ip) :: j  !! counter in \(y\)
    integer(ip) :: k  !! counter in \(z\)
    integer(ip) :: q  !! counter in \(q\)
    integer(ip) :: r  !! counter in \(r\)
    integer(ip) :: jf  !! output index

    f = 0.0_wp

    if (size(f,kind=ip)/=size(bcoef,6,kind=ip)) then
        iflag = 6002_ip
        return
    end if
    if (size(bcoef,1,kind=ip)/=nx .or. &
        size(bcoef,2,kind=ip)/=ny .or. &
        size(bcoef,3,kind=ip)/=nz .or. &
        size(bcoef,4,kind=ip)/=nq .or. &
        size(bcoef,5,kind=ip)/=nr) then
        iflag = 6004_ip
        return
    end if

    iflag = check_value(xval,tx,1_ip,extrap); if (iflag/=0_ip) return
    iflag = check_value(yval,ty,2_ip,extrap); if (iflag/=0_ip) return
    iflag = check_value(zval,tz,3_ip,extrap); if (iflag/=0_ip) return
    iflag = check_value(qval,tq,4_ip,extrap); if (iflag/=0_ip) return
    iflag = check_value(rval,tr,5_ip,extrap); if (iflag/=0_ip) return

    ! interval and basis function weights in each direction (once for all outputs)
    call dbvalu_weights(tx,nx,kx,idx,xval,inbvx,leftx,bx,iflag,extrap)
    if (iflag/=0_ip) return
    call dbvalu_weights(ty,ny,ky,idy,yval,iloy,lefty,by,iflag,extrap)
    if (iflag/=0_ip) return
    call dbvalu_weights(tz,nz,kz,idz,zval,iloz,leftz,bz,iflag,extrap)
    if (iflag/=0_ip) return
    call dbvalu_weights(tq,nq,kq,idq,qval,iloq,leftq,bq,iflag,extrap)
    if (iflag/=0_ip) return
    call dbvalu_weights(tr,nr,kr,idr,rval,ilor,leftr,br,iflag,extrap)
    if (iflag/=0_ip) return

    ix = leftx - kx
    iy = lefty - ky
    iz = leftz - kz
    iq = leftq - kq
    ir = leftr - kr

    do jf = 1_ip, size(f,kind=ip)
        ! contract the coefficients one direction at a time
        do r=1_ip,kr
            do q=1_ip,kq
                do k=1_ip,kz
                    do j=1_ip,ky
                        w4(j,k,q,r) = dot_product(bx,bcoef(ix+1_ip:leftx,iy+j,iz+k,iq+q,ir+r,jf))
                    end do
                    w3(k,q,r) = dot_product(by,w4(:,k,q,r))
                end do
                w2(q,r) = dot_product(bz,w3(:,q,r))
            end do
            w1(r) = dot_product(bq,w2(:,r))
        end do
        f(jf) = dot_product(br,w1)
    end do

    end subroutine db5val_multi