initialize_6d_auto_knots Subroutine

    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))
    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,iflag)

    if (iflag==0_ip) then
        call me%set_extrap_flag(extrap)
    end if

    me%initialized = iflag==0_ip
    me%iflag = iflag

    end subroutine initialize_6d_auto_knots