Tag: module

Here is a (reduced version of) a module for Bsplines, which I have used multiple times for my work on atomic and molecular physics. I wrote it more than a decade ago, with the occasional tweak, and I think it would benefit from a review. It is a bit long, but it is quite useful: it can compute arbitrary matrix elements between B-splines of the form
$$ int dx (D^n B_i(x)) f(x;mathbf{Y}) (D^m B_j(x)) $$
where $B_i(x)$ and $B_j(x)$ are Bspline functions, $D^n$ represents the n-th derivative, and $f(x;mathbf{Y})$ is an arbitrary real function of real variable that takes an optional array of parameters $mathbf{Y}$. I illustrate some of the capabilities of the module with a short test at the end that computes the eigenvalues of a quantum harmonic oscillator.
The code can arguably be improved here and there, but otherwise it works fine and I hope it can be interesting to some.

The main question that I have is on how to pass the function in the integral. At the moment, I pass it as a pointer to a pure function with a specified interface, and I have to explicitly define that function in the code, which is extremely annoying. I would rather like to be able to define arbitrary anonymous functions (lambdas) at run time, but I am not aware of how this could be done in fortran. Any suggestion?

Note: the code require blas and lapack to be linked. If you have ifort and mkl, you can compile it with the options -O2 -mkl

!! ModuleBSplines, Copyright (C) 2020 Luca Argenti, PhD - Some Rights Reserved
!! ModuleBSplines is licensed under a
!! Creative Commons Attribution-ShareAlike 4.0 International License.
!! A copy of the license is available at <http://creativecommons.org/licenses/by-nd/4.0/>.
!!
!! Luca Argenti is Associate Professor of Physics, Optics and Photonics
!! Department of Physics and the College of Optics
!! University of Central Florida
!! 4111 Libra Drive, Orlando, Florida 32816, USA
!! email: luca.argenti@ucf.edu

module ModuleBSpline

  use, intrinsic :: ISO_FORTRAN_ENV

  implicit none

  private

  interface 
     Pure DoublePrecision function D2DFun( x, parvec ) 
       DoublePrecision          , intent(in) :: x
       DoublePrecision, optional, intent(in) :: parvec(*)
     end function D2DFun
  end interface
  public D2DFun

  ! {{{ Private Attributes

  !
  logical                    :: INITIALIZED = .FALSE.
  !
  integer        , parameter :: IOMSG_LENGTH= 100
  !.. Local Gaussian weights
  DoublePrecision, parameter :: PI = 3.14159265358979323846d0
  DoublePrecision, parameter :: GAUSS_POINTS_CONVERGENCE_THRESHOLD = 2.d-15
  integer        , parameter :: NGAUSS = 64
  DoublePrecision            :: Gauss_points(NGAUSS)
  DoublePrecision            :: Gauss_weight(NGAUSS)
  !
  !.. Local Factorial Parameters and variables
  integer, parameter :: DIMFACT = 127
  DoublePrecision    :: fact(0:DIMFACT)
  DoublePrecision    :: tar(0:((DIMFACT+1)*(DIMFACT+2)/2-1))
  DoublePrecision    :: parfactMat(0:DIMFACT,0:DIMFACT)
  !
  !.. B-spline parameters
  integer        , parameter :: MAX_NUMBER_NODES = 10000
  DoublePrecision, parameter :: NOD_THRESHOLD    = 1.d-15
  integer        , parameter :: MAXNR            = 20
  !

  ! }}}

  !> Set of B-splines (spline basis functions)
  !!
  ! {{{ Detailed Description

  !> The spline are defined as piecewise polynomials of order f$kf$ 
  !! (maximum degree f$k-1f$), f$mathcal{C}^inftyf$ everywhere, except at a fixed set of knots f${t_i}f$ 
  !! where they are at least f$mathcal{C}^{k_i-2}f$, with f$k geq k_igeq 1f$. 
  !! f$nu_i=k-k_i+1f$ is called the knot multiplicity because the same spline basis can be obtained 
  !! from a basis with f$sum_i nu_if$ maximally regular knots in the limit where sets of f$nu_if$ 
  !! contiguous knots coalesce to f$if$-th knots. In other terms, the spline space is specified by an order 
  !! f$kf$ and a non-decreasing set of knots. The total dimension of spline space is f$n+kf$, but if 
  !! we specialize to the subset which is zero below the lowest and above the highest knots we are left 
  !! with a space of dimension f$n-kf$. As a consequence every spline extends at least over f$kf$ adjacent 
  !! intervals (f$k+1f$ consecutive knots, when counted with their multiplicity). The f$n-kf$ independent 
  !! splines restricted to f$kf$ intervals, clearly a basis of the spline space, are commonly called B-splines.
  !! Following deBoorcite{deBoor}, we define the f$if$-th B-spline f$B_i^k(x)f$ of order f$kf$ and which 
  !! extends from the knot f$t_if$ to the knot f$t_{i+k}f$, as follows:
  !! f{eqnarray}
  !! B_i^1(x)&=&theta(x-t_i),cdot,theta(t_{i+1}-x)\
  !! B_i^k(x)&=&frac{x-t_i}{t_{i+k-1}-t_i},,
  !!            B_i^{k-1}(x)+frac{t_{i+k}-x}{t_{i+k}-t_{i+1}},,
  !!            B_{i+1}^{k-1}(x).
  !! f}
  !! In the following, unless otherwise stated, we shall refer to standard set of knots where the 
  !! first and last breakpoints are f$kf$ times degenerate, while the other nodes are non degenerate:
  !! f(
  !! t_1=t_2=ldots=t_kleq t_{k+1}leqldotsleq t_{n-k}leq t_{n-k+1}=t_{n-k+2}=ldots=t_{n}
  !! f)
  !! The use of B-splines in atomic and molecular physics calculations is reviewed in 
  !! (cite{rpp.64.1815})(http://iopscience.iop.org/0034-4885/64/12/205/). 
  !! 
  !! B-splines are invariant upon affine transformations: that is, if
  !! f{equation}
  !! mathrm{if}quad xto x'=a,x +b,quad t_ito t_i'=a,t_i+bquadmathrm{then}
  !! quad{B_i^{k}}'(x')=B_i^k(x).
  !! f}
  !! 
  !! It is useful to define the f$L^2f$-normalized B-splines as
  !! f{equation}
  !! bar{B}_i(x)equiv B_i(x)/|B_i|_{L^2}
  !! f}
  !! 
  ! }}} 
  type, public :: ClassBSpline
     !
     private
     !> Number of non-degenerate nodes.
     integer :: NNodes
     !
     !> Spline order.
     !! A B-spline of order f$ kf$ is a piecewise polynomials
     !! of order f$k-1f$ which extends over f$kf$ intervals (including
     !! degenerate intervals, i.e., those with zero length).
     !! The quality of the Bspline evaluation is severely compromised
     !! for high order. At order around 35, in the case of just two nodes
     !! (which is the most demanding condition), the numerical error starts
     !! being of the same order of magnitude of the function itself.
     !! For order f$simf$20, there are occasional points of irregularity
     !! for a reasonable number of intervals (e.g.,f$geq 10f$). Anyhow,
     !! these caveats arise only in cases which already lie outside
     !! the range of applicability of the BSplines as a variational basis, 
     !! which is above order 12 or so.
     integer :: Order
     !
     !> Total number of B-spline. 
     !! f(
     !! NBsplines = NNodes + Order - 2
     !!f)
     integer :: NBSplines
     !
     !> grid which contains the node positions. 
     !! - f$g_{-k+1},=,g_{-k+2},=,cdots,=,g_{-1},=,g_0f$ 
     !!   is the first f$ kf$ times degenerate node.
     !! - f$g_{n-1},=,g_{n},=,cdots,=,g_{n+k-2}f$ is the 
     !!   last f$ kf$ times degenerate node.
     !! The set of distinct nodes is mapped to the positions (0:n-1) of grid;  
     !! and hence, the upper bound of the i-th BSpline support is grid(i)
     real(kind(1d0)),  allocatable :: grid(:)
     !
     !> Normalization constants.
     !! f(
     !!    mathrm{f}_i=|B_i|_{L^2}^{-1}=frac{1}{ sqrt{ langle B_i | B_i rangle } }
     !! f)
     real(kind(1d0)), private, allocatable :: f(:)
     !
     !
     !> Matrix of polynomial coefficients.
     !! The support $mathcal{D}_i$ of the f$if$-th B-spline comprises f$kf$ consecutive intervals:
     !! f(
     !!      mathcal{D}_i = (g_{i-k},g_i) = (g_{i-k+1},g_{i-k+1}) cup cdots cup (g_{i-1}:g_i).
     !! f)
     !! If we enumerate the individual intervals from f$0f$ to f$k-1f$, the explicit expression of 
     !! the f$if$-th B-spline in interval $j$ is
     !! f(
     !!     B_i(x) = sum_{m=0}^{k-1}
     !!        left(
     !!                frac{x-g_{i-k+j}}{g_{i-k+j+1}-g_{i-k+j}}
     !!        right)^m ,,  c(m,j,i)
     !! f)
     real(kind(1d0)), private, allocatable :: c(:,:,:)
     !
     !> Indicates whether ClassBSpline has been initialized or not.
     logical, private :: INITIALIZED=.FALSE.
     !
   contains

     !.. Basics
     generic   :: Init =>  BSplineInit
     procedure :: Free =>  BSplineFree
     procedure :: Save =>  ClassBSpline_Save
     procedure :: Load =>  ClassBSpline_Load

     !.. Accessors
     !> Gets the number of nodes.
     procedure :: GetNNodes         =>  BSplineGetNNodes
     !> Gets the B-splines order.
     procedure :: GetOrder          =>  BSplineGetOrder
     !> Gets the number of B-splines.
     procedure :: GetNBSplines      =>  BSplineGetNBSplines
     !> Gets the position of a requested node.
     procedure :: GetNodePosition   =>  BSplineNodePosition
     !> Gets the normalization factor of a requested B-spline.
     procedure :: GetNormFactor     =>  BSplineNormalizationFactor

     !> Evaluates either a choosen B-spline function, 
     !! or a linear combination of B-splines, both evaluated in certain position. 
     generic   :: Eval              =>  BSplineEval, BSplineFunctionEval

     !> Tabulates in a slected domain either  a choosen B-spline function, 
     !! or a linear combination of B-splines.
     generic   :: Tabulate          =>  BSplineTabulate, BSplineFunctionTabulate

     !> Computes the integral
     !! f(
     !!    int_{a}^{b} r^{2}dr
     !!       frac{d^{n_{1}}Bs_{i}(r)}{dr^{n_{1}}} f(r)
     !!       frac{d^{n_{2}}Bs_{j}(r)}{dr^{n_{2}}}
     !!f)
     !! Where f$f(r)f$ is a local operator, and if a break point f$BPf$ is introduced, 
     !! then the integral is splitted in two parts
     !! f(
     !!    int_{a}^{BP} r^{2}dr
     !!       frac{d^{n_{1}}Bs_{i}(r)}{dr^{n_{1}}} f(r)
     !!       frac{d^{n_{2}}Bs_{j}(r)}{dr^{n_{2}}} + 
     !!    int_{BP}^{b} r^{2}dr
     !!       frac{d^{n_{1}}Bs_{i}(r)}{dr^{n_{1}}} f(r)
     !!       frac{d^{n_{2}}Bs_{j}(r)}{dr^{n_{2}}}
     !! f)
     generic   :: Integral =>  BSplineIntegral
     procedure :: BSplineIntegral

     !.. Assignment
     !> Copies the B-spline set from one ClassBSpline to another. 
     generic, public :: ASSIGNMENT(=) => CopyBSplineSet

     !.. Private Interface
     procedure, private :: BSplineInit
     procedure, private :: CopyBSplineSet
     procedure, private :: BSplineEval
     procedure, private :: BSplineFunctionEval
     procedure, private :: BSplineTabulate
     procedure, private :: BSplineFunctionTabulate

  end type ClassBSpline

contains

  subroutine ClassBspline_Save( self, uid )
    class(ClassBspline), intent(in) :: self
    integer            , intent(in)  :: uid
    integer :: i
    write(uid,"(*(x,i0))") &
         self%Nnodes, self%Order, self%NBsplines, &
         size(self%grid), size(self%f), &
         (size(self%c,i),i=1,3) 
    write(uid,"(*(x,e24.16))") self%grid
    write(uid,"(*(x,e24.16))") self%f
    write(uid,"(*(x,e24.16))") self%c
  end subroutine ClassBspline_Save
   
  subroutine ClassBspline_Load( self, uid )
    class(ClassBspline), intent(out) :: self
    integer            , intent(in)  :: uid
    integer :: ng, nf, nc1, nc2, nc3
    call self%free()
    read(uid,*) self%Nnodes, self%Order, self%NBsplines, &
         ng, nf, nc1, nc2, nc3
    allocate(self%grid(ng))
    allocate(self%f   (nf))
    allocate(self%c(nc1,nc2,nc3))
    read(uid,*) self%grid
    read(uid,*) self%f
    read(uid,*) self%c
  end subroutine ClassBspline_Load


  !> Initialize module basic constants
  subroutine Init_Spline_Module()
    call InitGaussPoints
    call InitFactorials
    INITIALIZED=.TRUE.
  end subroutine Init_Spline_Module


  subroutine CheckInitialization(s)
    Class(ClassBSpline), intent(in) :: s
    if(.not.s%INITIALIZED) error stop
  end subroutine CheckInitialization


  !> Initialize the points and weights for 
  !! Gauss integration
  subroutine InitGaussPoints()
    integer         :: i,j
    DoublePrecision :: p1, p2, p3, pp, z, z1
    do i=1,floor((NGAUSS+1)/2.d0)
       z=cos(PI*(i-.25d0)/(NGAUSS+.5d0))
       inna : do
          p1=1.d0
          p2=0.d0
          do j=1,NGAUSS
             p3=p2
             p2=p1
             p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j
          enddo
          pp=NGAUSS*(z*p1-p2)/(z*z-1.d0)
          z1=z
          z=z1-p1/pp
          if(abs(z-z1)<=GAUSS_POINTS_CONVERGENCE_THRESHOLD)exit inna
       enddo inna
       Gauss_Points(i)=(1.d0-z)/2.d0
       Gauss_Points(NGAUSS+1-i)=(1.d0+z)/2.d0
       Gauss_Weight(i)=1.d0/((1.d0-z*z)*pp*pp)
       Gauss_Weight(NGAUSS+1-i)=1.d0/((1.d0-z*z)*pp*pp)
    enddo
  end subroutine InitGaussPoints


  subroutine InitFactorials()
    integer :: i,j,k,n
    !
    !Initialize Factorials
    fact(0)=1.d0
    do i = 1, DIMFACT
       fact(i)=fact(i-1)*dble(i)
    enddo
    !
    !.. Initialize Binomials
    k=0
    tar=0.d0
    do i = 0,DIMFACT
       do j = 0,i
          tar(k)=fact(i)/fact(i-j)/fact(j)
          k=k+1
       enddo
    enddo
    !
    parfactmat=1.d0
    do n=1,DIMFACT
       do k=1,n
          do i=n-k+1,n
             parfactmat(n,k)=parfactmat(n,k)*dble(i)
          enddo
       enddo
    enddo
    !
  end subroutine InitFactorials

  integer function BSplineGetOrder(s) result(Order)
    class(ClassBSpline), intent(in) :: s
    Order=s%Order
  end function BSplineGetOrder

  integer function BSplineGetNNodes(s) result(NNodes)
    class(ClassBSpline), intent(in) :: s
    NNodes=s%NNodes
  end function BSplineGetNNodes

  integer function BSplineGetNBSplines(s) result(NBSplines)
    class(ClassBSpline), intent(in) :: s
    NBSplines=s%NBSplines
  end function BSplineGetNBSplines

  !> Return the position of a node
  DoublePrecision function BSplineNodePosition(SplineSet,Node) result(Position)
    Class(ClassBSpline), intent(in) :: SplineSet
    integer           , intent(in) :: Node
    if( Node <= 1 )then
       Position = SplineSet%Grid(0)
    elseif( Node >= SplineSet%NNodes )then
       Position = SplineSet%Grid( SplineSet%NNodes-1 )
    else
       Position = SplineSet%Grid( Node-1 )
    endif
    return
  end function BSplineNodePosition

  !> Return the position of a node
  DoublePrecision function BSplineNormalizationFactor(SplineSet,Bs) &
       result(Normalization)
    Class(ClassBSpline), intent(in) :: SplineSet
    integer            , intent(in) :: Bs
    Normalization = 0.d0
    if(Bs<1.or.Bs>SplineSet%NBSplines)return
    Normalization=SplineSet%f(Bs)
  end function BSplineNormalizationFactor


  !> Copies the B-spline set sOrigin to sDestination 
  subroutine CopyBSplineSet(sDestination,sOrigin)
    Class(ClassBSpline), intent(inout):: sDestination
    Class(ClassBSpline), intent(in)   :: sOrigin
    !
    sDestination%NNodes=sOrigin%NNodes
    sDestination%Order=sOrigin%Order
    sDestination%NBSplines=sOrigin%NBSplines
    !
    if(allocated(sDestination%Grid))deallocate(sDestination%Grid)
    allocate(sDestination%Grid(1-sDestination%Order:sDestination%NBSplines))
    sDestination%Grid=sOrigin%Grid
    !
    if(allocated(sDestination%c))deallocate(sDestination%c)
    allocate(sDestination%c(&
         0:sDestination%Order-1 ,&
         0:sDestination%Order-1 ,&
         1:sDestination%NBSplines+sDestination%Order-2))
    sDestination%c=sOrigin%c
    !
    if(allocated(sDestination%f))deallocate(sDestination%f)
    allocate(sDestination%f(sDestination%NBSplines))
    sDestination%f=sOrigin%f
    !
    sDestination%INITIALIZED=.TRUE.
  end subroutine CopyBSplineSet


  subroutine CheckNodesAndOrder(NumberOfNodes,Order,STAT)
    integer, intent(in)  :: NumberOfNodes, Order
    integer, intent(out) :: STAT
    STAT = -1; if ( NumberOfNodes < 2 ) return
    STAT = -2; if (    order      < 1 ) return
    STAT = 0
  end subroutine CheckNodesAndOrder
  

  !> Initialize a BSpline set
  Subroutine BSplineInit(s,NumberOfNodes,order,grid,IOSTAT)
    !
    !> Bspline basis set to be initialized
    Class(ClassBSpline), intent(inout)           :: s
    !> Number of Nodes
    integer            , intent(in)              :: NumberOfNodes
    !> Spline order
    integer            , intent(in)              :: order
    !> Grid of nodes (without repetition of degenerate nodes)
    DoublePrecision    , intent(in)              :: grid(1:NumberOfNodes)
    !> IOSTAT = 0 on successful exit, IOSTAT/=0 otherwise
    !! If IOSTAT is not specified, the programs terminates with
    !! an error.
    integer            , intent(out),   optional :: IOSTAT

    !Local variables
    integer                    :: i, Bs, Status
    procedure(D2DFUN), pointer :: FunPtr

    if(present(IOSTAT)) IOSTAT=0
    if(.not.INITIALIZED)call Init_Spline_Module()

    !.. Check input
    call CheckNodesAndOrder(NumberOfNodes,Order,Status)
    if(Status/=0)then
       if(present(IOSTAT))then
          IOSTAT=Status
          return
       endif
       error stop 
    endif

    s%NNodes    = NumberOfNodes
    s%Order     = order
    s%NBSplines = NumberOfNodes + order - 2

    !.. Check Grid format
    if(  LBOUND(Grid,1) > 1             .or. &
         UBOUND(Grid,1) < NumberOfNodes ) then
       if(present(IOSTAT))then
          IOSTAT=-4
          return
       endif
       error stop
    endif
    do i=2,NumberOfNodes
       if( Grid(i) < Grid(i-1) )then
          if(present(IOSTAT))then
             IOSTAT=1
             return
          else
             error stop
          endif
       endif
    enddo

    !.. Maps the Grid to (0:n-1); hence, the upper bound
    !   of the i-th BSpline support is g(i)
    if( allocated( s%Grid ) ) deallocate( s%Grid )
    allocate( s%Grid( -(order-1) : NumberOfNodes-1 + order-1 ) )
    s%Grid( 1-order : -1 ) = Grid(1)
    s%Grid( 0:NumberOfNodes-1 ) = Grid(1:NumberOfNodes)
    s%Grid( NumberOfNodes:NumberOfNodes+order-2 ) = Grid(NumberOfNodes)

    !.. Initializes the BSpline coefficients normalized
    !   so that the B-spline set is a partition of unity
    !..
    if( allocated( s%c ) ) deallocate( s%c )
    allocate( s%c( 0:order-1, 0:order-1, 1:s%NBSplines ) )
    s%c = 0.d0
    do Bs=1,s%NBSplines
       call ComputeCoefficientsSingleBSpline(s%Order,s%Grid(Bs-s%Order),s%c(0,0,Bs))
    enddo
    s%INITIALIZED=.TRUE.

    !.. Initialize the normalization factors
    !   ||Bs_i(x)*s%f(i)||_L2=1.d0
    !..
    if( allocated( s%f ) ) deallocate( s%f )
    allocate( s%f( 1 : s%NBSplines ) )
    s%f = 0.d0
    funPtr=>Unity
    do i = 1, s%NBSplines
       s%f(i)=1.d0/sqrt(s%Integral(FunPtr,i,i))
    enddo

  end Subroutine BSplineInit


  Pure DoublePrecision function Unity(x,parvec) result(y)
    DoublePrecision, intent(in) :: x
    DoublePrecision, optional, intent(in) :: parvec(*)
    y=1.d0
  end function Unity


  !> Given the order and a set of order+1 consecutive nodes,
  !> builds the corresponding (uniquely defined) B-spline.
  subroutine ComputeCoefficientsSingleBSpline(order,vec,Coefficients)
    !
    !> Order of the B-spline = number of (possibly degenerate) 
    !> intervals that form the B-spline support = Polynomials degree + 1
    !> (e.g., order=4 means piecewise cubic polynomials, which require
    !> four parameters: f$ a_0 + a_1 x + a_2 * x^2 + a_3 * x^3 f$.
    integer        , intent(in)  :: order
    !
    !> vector of order+1 nodes. For Order = 5, for example, vec contains 
    !> the six bounds of the five consecutive intervals that form the 
    !> support of the B-spline, indexed from ONE to SIX.
    !
    !  vec(1)  vec(2)    vec(3)  vec(4) vec(5)    vec(6)
    !    |-------|---------|-------|------|---------|
    DoublePrecision, intent(in)  :: vec(1:order+1)
    !
    !> Matrix of polynomial coefficients in each interval.
    DoublePrecision, intent(out) :: Coefficients(0:order-1,0:order-1)
    !
    !Local variables
    integer         :: k, inter, tempBs
    DoublePrecision :: SpanBSpline
    DoublePrecision :: SpanCurrentInterval
    DoublePrecision :: SpanPriorToInterval
    DoublePrecision :: SpanToTheEnd
    DoublePrecision :: FactorConstant
    DoublePrecision :: FactorMonomial
    DoublePrecision, allocatable :: coe(:,:,:,:)
    !
    !.. Computes the B-spline in each available
    !   sequence of order consecutive intervals
    !
    !
    allocate( coe(0:order-1,0:order-1,order,order) )
    coe=0.d0
    coe(0,0,1,:)=1.d0
    !
    !.. Cycle over the spline order
    !
    do k = 2, order
       !
       !.. Cycle over the (temporary) Bsplines available 
       !   for any given order, which must all be computed
       !
       do tempBs = 1, order - k + 1
          !
          !.. Each B-spline of order k is obtained as the sum of 
          !   the contribution from the two B-splines of order k-1
          !   whose support is comprised in that of the final 
          !   B-spline.
          !
          !.. Contribution of the first B-spline of order k-1
          !   counted only if the B-spline is not degenerate
          !..
          SpanBSpline = vec(tempBs+k-1) - vec(tempBs)
          if( SpanBSpline > NOD_THRESHOLD )then
             do inter=0,k-2
                !
                !.. Constant Component
                SpanPriorToInterval = vec(tempBs+inter) - vec(tempBs)
                FactorConstant = SpanPriorToInterval / SpanBSpline 
                coe(0:k-2,inter,k,tempBs) = coe(0:k-2,inter,k,tempBs) + &
                     coe(0:k-2,inter,k-1,tempBs) * FactorConstant
                !
                !.. Monomial Component
                SpanCurrentInterval = vec(tempBs+inter+1) - vec(tempBs+inter)
                FactorMonomial = SpanCurrentInterval / SpanBSpline
                coe(1:k-1,inter,k,tempBs) = coe(1:k-1,inter,k,tempBs) + &
                     coe(0:k-2,inter,k-1,tempBs) * FactorMonomial
                !
             enddo
          endif
          !
          !.. Contribution of the second B-spline of order k-1
          !   counted only if the B-spline is not degenerate
          !..
          SpanBSpline = vec(tempBs+k) - vec(tempBs+1)
          if( SpanBSpline > NOD_THRESHOLD )then
             do inter=1,k-1
                !
                !.. Constant Component
                SpanToTheEnd = vec(tempBs+k)-vec(tempBs+inter)
                FactorConstant = SpanToTheEnd / SpanBSpline
                coe(0:k-2,inter,k,tempBs) = coe(0:k-2,inter,k,tempBs) + &
                     coe(0:k-2,inter-1,k-1,tempBs+1) * FactorConstant
                !
                !.. Monomial Component
                SpanCurrentInterval = vec(tempBs+inter+1) - vec(tempBs+inter)
                FactorMonomial = SpanCurrentInterval / SpanBSpline
                coe(1:k-1,inter,k,tempBs) = coe(1:k-1,inter,k,tempBs) - &
                     coe(0:k-2,inter-1,k-1,tempBs+1) * FactorMonomial
                !
             enddo
          endif
          !
       enddo
    enddo
    !
    Coefficients=coe(:,:,order,1)
    deallocate(coe)
    !
  end subroutine ComputeCoefficientsSingleBSpline


  !> Deallocates a ClassBSpline variable
  subroutine BSplineFree(s)
    Class(ClassBSpline), intent(inout) :: s
    if(allocated(s%Grid))deallocate(s%Grid)
    if(allocated(s%f))deallocate(s%f)
    if(allocated(s%c))deallocate(s%c)
    s%NNodes=0
    s%Order=0
    s%NBSplines=0
    s%INITIALIZED=.FALSE.
    return
  end subroutine BSplineFree


  !> Returns n if x is in f$(g_n,g_{n+1})f$.
  integer function which_interval(x,s)
    !
    DoublePrecision   , intent(in) :: x
    type(ClassBSpline), intent(in) :: s
    !
    integer, save :: i1 = 0
    integer       :: i2
    !
    if(x>s%Grid(i1).and.x<=s%Grid(i1+1))then
       which_interval=i1
       return
    endif
    if(x>s%Grid(i1+1).and.x<=s%Grid(min(i1+2,s%NNodes+s%Order-1)))then
       which_interval=i1+1
       return
    endif
    i1=0
    i2=0
    !
    if(x<s%Grid(0))then
       which_interval=-1
       return
    elseif(x>s%Grid(s%NNodes-1))then
       which_interval=s%NNodes!-1
       return
    endif
    if(x<s%Grid(i1))i1=0
    if(x>s%Grid(i2))i2=s%NNodes-1
    do while(i2-i1>1)
       if(x>s%Grid((i2+i1)/2))then
          i1=(i1+i2)/2
       else
          i2=(i1+i2)/2
       endif
    enddo
    which_interval=i1
    return
  end function which_interval


  !> Computes the n-th derivative of B-spline Bs, f$ n=0,1,2,ldotsf$
  !!
  ! {{{ Detailed Description:

  !! ---------------------
  !! B-splines are positive: expansion coefficients of represented functions are 
  !! akin to the functions themselves. Without rapid oscillations of coefficients, 
  !! the cancellation errors are kept at minimum. 
  !! In particular cite{deBoor},
  !! if f$f=sum_i B_ic_if$
  !! f{equation}
  !! left|c_i-frac{M+m}{2}right|leq D_{k}frac{M-m}{2},quad m
  !!  =min_{xin(a,b)} f(x),,,M=max_{xin(a,b)}f(x)
  !! f}
  !! where f$D_kf$ is a constant that depends only on f$kf$ and not on the
  !! particular partition of the f$(a,b)f$ interval.
  !! At any point, the evaluation of a linear combination of B-splines 
  !! requires the evaluation of f$kf$ basis functions only. 

  ! }}}
  DoublePrecision function BSplineEval(s,x,Bs,n_)
    Class(ClassBSpline), intent(in) :: s
    DoublePrecision    , intent(in) :: x
    integer            , intent(in) :: Bs
    integer, optional  , intent(in) :: n_
    !
    integer         :: i,j,k,n
    DoublePrecision :: r,a,w,OneOverInterval,OneOverIntervalToTheN
    !
    call CheckInitialization(s)
    !
    BSplineEval=0.d0
    if(Bs<1.or.Bs>s%NBSplines)return
    !
    n=0; if(present(n_)) n=n_
    if(n<0) error stop 
    if(n>=s%Order)return
    !
    i=which_interval(x,s)
    if(  i<max(0,Bs-s%Order).or.&
         i>min(s%NNodes-1,Bs-1))return
    !
    j=i-Bs+s%Order
    OneOverInterval=1.d0/(s%Grid(i+1)-s%Grid(i))
    r=(x-s%Grid(i))*OneOverInterval
    !
    a=1.d0
    w=0.d0
    do k=n,s%Order-1
       w=w+a*parfactmat(k,n)*s%c(k,j,Bs)
       a=a*r
    enddo
    !
    OneOverIntervalToTheN=1.d0
    do k=1,n
       OneOverIntervalToTheN=OneOverIntervalToTheN*OneOverInterval
    enddo
    BSplineEval=w*OneOverIntervalToTheN
    !
    return
  end function BSplineEval

  !> Computes the k-th derivative of a function expressed in terms of B-splines.
  DoublePrecision function BSplineFunctionEval(s,x,fc,k,SKIP_FIRST,Bsmin,Bsmax)
    class(ClassBSpline), intent(in) :: s
    DoublePrecision    , intent(in) :: x
    DoublePrecision    , intent(in) :: fc(:)
    integer, optional  , intent(in) :: k
    logical, optional  , intent(in) :: SKIP_FIRST
    integer, optional  , intent(in) :: Bsmin,Bsmax
    integer :: Bsmi, Bsma
    integer :: Bs,i,nskip
    call CheckInitialization(s)
    BSplineFunctionEval=0.d0
    i=which_interval(x,s)
    if(i<0.or.i>s%NNodes-2)return
    Bsmi=1
    nskip=0
    if(present(SKIP_FIRST))then
       if(SKIP_FIRST)then
          Bsmi=2
          nskip=1
       endif
    endif
    if(present(Bsmin))then
       Bsmi = max(Bsmin,Bsmi)
       nskip= max(Bsmi-1,nskip)
    endif
    Bsmi=max(Bsmi,i+1)
    Bsma = min(s%NBSplines,i+s%Order)
    if(present(Bsmax)) Bsma = min(Bsmax,Bsma)
    do Bs=Bsmi,Bsma
       BSplineFunctionEval=BSplineFunctionEval+s%Eval(x,Bs,k)*fc(Bs-nskip)
    enddo
    return
  end function BSplineFunctionEval
  
  !> Computes the  k-th derivatives of a choosen B-spline evaluated in a position array. 
  subroutine BSplineTabulate(s,ndata,xvec,yvec,Bs,n_) 
    Class(ClassBSpline), intent(in) :: s
    integer            , intent(in) :: ndata
    DoublePrecision    , intent(in) :: xvec(1:ndata)
    DoublePrecision    , intent(out):: yvec(1:ndata)
    integer            , intent(in) :: Bs
    integer, optional  , intent(in) :: n_
    !
    integer :: i,n
    n=0;if(present(n_))n=n_
    do i=1,ndata
       yvec(i)=s%Eval(xvec(i),Bs,n)
    enddo
  end subroutine BSplineTabulate

  !> Computes the  k-th derivatives of a function expressed in terms of B-splines and evaluated in an array of positions. 
  subroutine BSplineFunctionTabulate(s,ndata,xvec,yvec,FunVec,n_,SKIP_FIRST) 
    Class(ClassBSpline), intent(in) :: s
    integer            , intent(in) :: ndata
    DoublePrecision    , intent(in) :: xvec(1:ndata)
    DoublePrecision    , intent(out):: yvec(1:ndata)
    DoublePrecision    , intent(in) :: FunVec(:)
    integer, optional  , intent(in) :: n_
    logical, optional  , intent(in) :: SKIP_FIRST
    logical :: SKIP_FIRST_LOC
    !
    integer :: i,n
    SKIP_FIRST_LOC=.FALSE.
    if(present(SKIP_FIRST))then
       SKIP_FIRST_LOC=SKIP_FIRST
    endif
    n=0;if(present(n_))n=n_
    do i=1,ndata
       yvec(i)=s%Eval(xvec(i),FunVec,n,SKIP_FIRST=SKIP_FIRST_LOC)
    enddo
  end subroutine BSplineFunctionTabulate


  !> Computes the integral
  !! f(
  !!    int_{a}^{b} frac{d^{n_{1}}Bs_{i}(r)}{dr^{n_{1}}} f(r)
  !!                 frac{d^{n_{2}}Bs_{j}(r)}{dr^{n_{2}}} r^{2}dr
  !! f)
  !! Where f$f(r)f$ is a local operator, and if a break point f$BPf$ 
  !! is introduced, then the integral is splitted in two parts
  !! f(
  !!    int_{a}^{BP} frac{d^{n_{1}}Bs_{i}(r)}{dr^{n_{1}}} f(r)
  !!                  frac{d^{n_{2}}Bs_{j}(r)}{dr^{n_{2}}} r^{2}dr + 
  !!    int_{BP}^{b} frac{d^{n_{1}}Bs_{i}(r)}{dr^{n_{1}}} f(r)
  !!                  frac{d^{n_{2}}Bs_{j}(r)}{dr^{n_{2}}}r^{2}dr
  !! f)
  DoublePrecision function BSplineIntegral( &
       s                 , &
       FunPtr            , & 
       Bs1               , &
       Bs2               , &
       BraDerivativeOrder, &
       KetDerivativeOrder, &
       LowerBound        , &
       UpperBound        , &
       BreakPoint        , &
       parvec            ) &
       result( Integral )
    !
    Class(ClassBSpline)      , intent(in) :: s
    procedure(D2DFun)        , pointer    :: FunPtr
    integer                  , intent(in) :: Bs1
    integer                  , intent(in) :: Bs2
    integer        , optional, intent(in) :: BraDerivativeOrder
    integer        , optional, intent(in) :: KetDerivativeOrder
    DoublePrecision, optional, intent(in) :: LowerBound
    DoublePrecision, optional, intent(in) :: UpperBound
    DoublePrecision, optional, intent(in) :: BreakPoint
    DoublePrecision, optional, intent(in) :: Parvec(*)
    !
    integer         :: n1, n2, iostat
    DoublePrecision :: a, b
    real(kind(1d0)) :: NewBreakPoint
    !
    Integral=0.d0
    call CheckInitialization(s)
    call CheckParameters( iostat )
    if( iostat/=0 )return
    if( present(BreakPoint) )then
       !***
       if ( BreakPoint < a ) then
          NewBreakPoint = a + epsilon(1d0)
       elseif (  BreakPoint > b ) then
          NewBreakPoint = b - epsilon(1d0)
       end if
       !***
       Integral = &
            BSplineDriverIntegral(s,FunPtr,Bs1,Bs2,n1,n2, a, BreakPoint,    Parvec ) + &
            BSplineDriverIntegral(s,FunPtr,Bs1,Bs2,n1,n2, NewBreakPoint, b, Parvec )
    else
       Integral = BSplineDriverIntegral(s,FunPtr,Bs1,Bs2,n1,n2,a,b, Parvec )
    endif
    !
    return
    !
  contains
    !
    subroutine CheckParameters( IOStat )
      integer, intent(out) :: IOStat
      call CheckBSplineIndexes( IOStat ); if( IOStat /= 0 ) return
      call CheckIntegralBounds( IOStat ); if( IOStat /= 0 ) return
      call CheckDerivativeOrder
    end subroutine CheckParameters
    !
    subroutine CheckBSplineIndexes( IOStat )
      integer, intent(out) :: IOStat
      IOStat=1
      if( min(Bs1,Bs2) <  1           ) return      
      if( max(Bs1,Bs2) >  s%NBSplines ) return
      if( abs(Bs1-Bs2) >= s%Order     ) return
      IOStat=0
    end subroutine CheckBSplineIndexes
    !
    subroutine CheckIntegralBounds( IOStat )
      integer, intent(out) :: IOStat
      IOStat=1
      a = s%Grid( max(Bs1,Bs2) - s%Order )
      b = s%Grid( min(Bs1,Bs2) )
      if(present(LowerBound))then
         if( LowerBound >= b )return
         a=max(a,LowerBound)
      endif
      if(present(UpperBound))then
         if( UpperBound <= a )return
         b=min(b,UpperBound)
      endif
      IOStat=0
    end subroutine CheckIntegralBounds
    !
    subroutine CheckDerivativeOrder
      n1=0; if(present(BraDerivativeOrder)) n1=BraDerivativeOrder
      n2=0; if(present(KetDerivativeOrder)) n2=KetDerivativeOrder
    end subroutine CheckDerivativeOrder
    !
  end function BSplineIntegral


  ! {{{ Detailed Description

  !> Compute the integral
  !! f(
  !!    int_a^b dr 
  !!    frac{d^{n_1}B_1(r)}{dr^{n_1}} f(r)
  !!    frac{d^{n_2}B_2(r)}{dr^{n_2}},
  !! f)
  !! where f$f(r)f$ is a *smooth* function in the
  !! intersection of the supports of the two B-splines.
  !! If the first (last) boundary of the integration 
  !! interval is not specified, then the value a=-oo
  !! (b=+oo) is assumed.
  !!
  !!
  !! Properties of B-spline integrals:
  !! ---------------------------------
  !! Since the support f$D_if$ of the f$if$-th B-spline is formed by
  !! f$kf$ consecutive intervals, the integrals between two B-splines and a 
  !! local operator are zero unless their indices differ less than the order f$kf$:
  !! f{equation}
  !! langle B_i|O|B_jrangle=int_{D_i cap D_j} B_i(x) o(x) B_j(x) dx.
  !! f}
  !! As a consequence, matrices are sparse and often band-diagonal. 
  !! On uniform grids, matrices of translationally invariant operators, that 
  !! is with kernel f$o(x,y)=o(x-y)f$, are Toeplitz:
  !! f{eqnarray}
  !! langle B_i|O|B_jrangle=langle B_{i+n}|O|B_{j+n}rangle.
  !! f}
  !! If an hermitian operator f$Of$ is both local and translationally invariant, 
  !! like the identity and the kinetic energy, its matrix on a uniform grid
  !! is Toeplitz and banded, in other terms it is defined by just f$kf$ numbers. 
  !! 
  !! With Gauss-type integration technique, the matrix elements of operators 
  !! with polynomial kernels are exact. The error f$epsilonf$ in the approximation
  !! of a f$mathcal{C}^kf$ function f$ff$ with B-splines, is bounded by
  !! f(
  !!     epsilonleq mathrm{const}_k|mathbf{t}|^k|D^kf|,quad |mathbf{t}|=
  !!         max_i(t_{i+1}-t_i),quad |f|=max_{xin(a,b)}|f(x)|
  !! f)
  !!

  ! }}}
  DoublePrecision function BSplineDriverIntegral(s,FunPtr,Bs1,Bs2,n1,n2,a,b,parvec) &
       result( Integral )
    !
    !.. Assumes that the input data have been sanitized
    !
    Class(ClassBSpline), intent(in) :: s
    procedure(D2DFun)  , pointer    :: FunPtr
    integer            , intent(in) :: Bs1, Bs2
    integer            , intent(in) :: n1, n2
    DoublePrecision    , intent(in) :: a, b
    DoublePrecision, optional, intent(in) :: parvec(*)
    !
    integer         :: Interval, IntervalMin, IntervalMax
    DoublePrecision :: LowerBound, UpperBound
    !
    DoublePrecision :: PartialIntegral
    DoublePrecision :: IntervalWidth
    DoublePrecision :: Radius
    integer         :: iGauss
    !
    DoublePrecision :: aPlus, bMinus
    !
    Integral=0.d0
    !
    aPlus  = UpperLimitTo( a )
    bMinus = LowerLimitTo( b )
    IntervalMin = which_interval( aPlus,  s )
    IntervalMax = which_interval( bMinus, s )
    !
    do Interval = IntervalMin, IntervalMax
       LowerBound = max( a, s%Grid( Interval   ) )
       UpperBound = min( b, s%Grid( Interval+1 ) )
       !
       IntervalWidth = UpperBound - LowerBound
       if( IntervalWidth < NOD_THRESHOLD ) cycle
       !
       PartialIntegral = 0.d0
       do iGauss=1,NGAUSS
          !
          Radius = LowerBound + IntervalWidth * Gauss_Points( iGauss )
          PartialIntegral = PartialIntegral + &
               FunPtr(Radius, Parvec) * &
               s%Eval(Radius,Bs1,n1)  * &
               s%Eval(Radius,Bs2,n2)  * &
               Gauss_weight( iGauss )
          !
       enddo
       PartialIntegral = PartialIntegral * IntervalWidth
       !
       Integral = Integral + PartialIntegral
       !
    enddo
    !
  end function BSplineDriverIntegral

  Pure DoublePrecision function UpperLimitTo(x) result(xPlus)
    DoublePrecision, intent(in) :: x
    DoublePrecision :: tin,eps
    eps=epsilon(1.d0)
    tin=tiny(x)*(1.d0+eps)
    xPlus=(x+tin)*(1.d0+eps)
  end function UpperLimitTo

  Pure DoublePrecision function LowerLimitTo(x) result(xMinus)
    DoublePrecision, intent(in) :: x
    DoublePrecision :: tin,eps
    eps=epsilon(1.d0)
    tin=tiny(x)*(1.d0+eps)
    xMinus=(x-tin)*(1.d0-eps)
  end function LowerLimitTo

end module ModuleBSpline


!> Computes the eigenvalues of a number of states of the quantum
!! harmonic oscillator, and compares them with the analytically
!! known result, $ E_n = n + 1/2 $ (when the parameters are chosen
!! such that hbar and the angular frequency are one.
Program TestModuleBSpline

  use, intrinsic :: iso_fortran_env, only : ERROR_UNIT, OUTPUT_UNIT
  use ModuleBSpline

  implicit none

  !.. Bspline parameters and variables
  integer        , parameter   :: BS_NNODS = 501
  integer        , parameter   :: BS_ORDER =  12
  real(kind(1d0)), parameter   :: BS_GRMIN = -20.d0
  real(kind(1d0)), parameter   :: BS_GRMAX =  20.d0
  real(kind(1d0)), parameter   :: BS_INTER =  BS_GRMAX - BS_GRMIN
  real(kind(1d0))              :: BS_GRID(BS_NNODS)
  type(ClassBSpline)           :: BSpline

  !.. Hamiltonian, overlap, and spectral-decomposition arrays
  real(kind(1d0)), allocatable :: Hmat(:,:), Smat(:,:)
  integer                      :: iNode, nEn, info

  !.. Initializes the BSpline set
  do iNode=1,BS_NNODS
     BS_GRID(iNode) = BS_GRMIN + &
          BS_INTER * dble(iNode-1) / dble(BS_NNODS-1)
  enddo
  call BSpline%Init( &
       BS_NNODS       , &
       BS_ORDER       , &
       BS_GRID        , &
       info           )
  if(info/=0) error stop

  !.. Due to the continuity requirement on the wavefunction,
  !   the first and last BSplines, which do not vanish at the
  !   boundary of the interval, cannot be used
  nEn = BSpline%GetNBsplines() - 2 

  BlockFillMatrices: block 
    integer         :: iBs1, iBs2
    real(kind(1d0)) :: Overlap, KineticEnergy, PotentialEnergy
    real(kind(1d0)) :: parvec(1)
    procedure(D2DFUN), pointer :: fPtrUni, fPtrPow
    fPtrUni => Unity
    fPtrPow => Power
    allocate(Smat(BS_ORDER,nEn))
    Smat=0.d0
    allocate(Hmat,source=Smat)
    parvec(1)=2.d0
    do iBs2 = 2, nEn + 1
       do iBs1 = max(2, iBs2-BS_ORDER+1), iBs2
          Overlap         = BSpline%Integral(fPtrUni,iBs1,iBs2)
          KineticEnergy   = BSpline%Integral(fPtrUni,iBs1,iBs2,1,1) / 2.d0
          PotentialEnergy = Bspline%Integral(fPtrPow,iBs1,iBs2,parvec=parvec) / 2.d0
          Smat(iBs1+BS_ORDER-iBs2,iBs2-1) = Overlap      
          Hmat(iBs1+BS_ORDER-iBs2,iBs2-1) = KineticEnergy + PotentialEnergy
       enddo
    enddo
  end block BlockFillMatrices

  BlockGeneralizedEigenproblem : Block
    real(kind(1d0)), allocatable :: Eval(:), Evec(:,:), Work(:)
    real(kind(1d0)), parameter   :: ERROR_THRESHOLD = 1.d-10 
    real(kind(1d0))              :: EigenvalueError
    integer                      :: iEn, info
    allocate(work(3*nEn),Eval(nEn),Evec(1,1))
    work=0.d0
    Eval=0.d0
    Evec=0.d0
    call DSBGV( 'N', 'U', nEn, BS_ORDER-1, BS_ORDER-1, &
         Hmat, BS_ORDER, Smat, BS_ORDER, Eval, Evec, 1, work, info )
    if( info /= 0 )then
       write(ERROR_UNIT,"(a,i0)") "DSBGV info : ", info 
       error stop
    endif
    deallocate(work)
    do iEn = 1, nEn
       EigenvalueError = Eval(iEn) - (iEn-0.5d0)
       if( EigenvalueError > ERROR_THRESHOLD )exit
       write(OUTPUT_UNIT,"(i4,x,e24.16)") iEn, EigenvalueError 
    enddo
    write(OUTPUT_UNIT,"(a,i0)") "Number of Accurate Eigenvalues : ",iEn-1
  end Block BlockGeneralizedEigenproblem
  
contains

  !.. Lambdas would be so useful ...
  Pure real(kind(1d0)) function Unity(x,parvec) result(y)
    DoublePrecision, intent(in) :: x
    DoublePrecision, optional, intent(in) :: parvec(*)
    y=1.d0
  end function Unity

  Pure real(kind(1d0)) function Power(x,parvec) result(res)
    real(kind(1d0)), intent(in) :: x
    real(kind(1d0)), optional, intent(in) :: parvec(*)
    res = x**parvec(1)
  end function Power
  
end Program TestModuleBSpline