NAG Library Manual, Mark 28.7
```!   G13NEF Example Program Text
!   Mark 28.7 Release. NAG Copyright 2022.
Module g13nefe_mod

!     G13NEF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Accessibility Statements ..
Private
Public                           :: chgpfn, get_data
Contains
Subroutine chgpfn(side,u,w,minss,v,cost,y,iuser,ruser,info)
!       Routine to calculate a proposed change point and associated cost
!       The cost is based on the likelihood of the gamma distribution

!       .. Use Statements ..
Use nag_library, Only: x07caf, x07cbf
!       .. Scalar Arguments ..
Integer, Intent (Inout)        :: info
Integer, Intent (In)           :: minss, side, u, w
Integer, Intent (Out)          :: v
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: cost(3)
Real (Kind=nag_wp), Intent (Inout) :: ruser(0:*), y(*)
Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
Real (Kind=nag_wp)             :: dn, shape, this_cost, tmp, ys
Integer                        :: floc, i, li, lloc
!       .. Local Arrays ..
Integer                        :: cexmode(3), texmode(3)
!       .. Intrinsic Procedures ..
Intrinsic                      :: log
!       .. Executable Statements ..

!       The gamma cost function used below can result in log(0) being taken
!       (if there is a segment of zeros in Y), this leads to a cost of -Inf
!       (which is correct), but we need to make sure that the compiler
!       doesn't stop at the creation of the -Inf

!       Save the current IEEE exception mode
Call x07caf(cexmode)

!       Set the IEEE exception mode to not trap division by zero
texmode(:) = cexmode(:)
texmode(2) = 0
Call x07cbf(texmode)

!       Extract shape from RUSER
shape = ruser(0)

!       Calculate the first and last positions for potential change
!       points, conditional on the fact that each sub-segment must be
!       at least MINSS wide
floc = u + minss - 1
lloc = w - minss

!       In order to calculate the cost of having a change point at I, we
!       need to calculate C(Y(FLOC:I)) and C(Y(I+1:LLOC)), where C(.) is
!       the cost function (based on the gamma distribution in this example).
!       Rather than calculate these values at each call to CHGPFN we store
!       the values for later use

!       If SIDE = 1 (i.e. we are working with a left hand sub-segment),
!       we already have C(Y(FLOC:I)) for this value of FLOC, so only need
!       to calculate C(Y(I+1:LLOC)), similarly when SIDE = 2 we only need
!       to calculate C(Y(FLOC:I))
!       When SIDE = -1, we need the cost of the full segment, which we do
!       in a forwards manner (calculating C(Y(FLOC:I)) in the process), so
!       when SIDE = 0 we only need to calculate C(Y(I:1:LLOC))

!       Get the intermediate costs
ys = 0.0_nag_wp
dn = 0.0_nag_wp
If (side==0 .Or. side==1) Then
!         RUSER(2*I) = C(Y(I+1:W))
Do i = w, floc + 1, -1
dn = dn + 1.0_nag_wp
tmp = dn*shape
ys = ys + y(i)
ruser(2*i-2) = 2.0_nag_wp*tmp*(log(ys)-log(tmp))
End Do

Else
!         RUSER(2*I-1) = C(Y(U:I))
If (side==-1) Then
li = w
Else
li = lloc
End If
Do i = u, li
dn = dn + 1.0_nag_wp
tmp = dn*shape
ys = ys + y(i)
ruser(2*i-1) = 2.0_nag_wp*tmp*(log(ys)-log(tmp))
End Do
End If

If (side>=0) Then
!         Need to find a potential change point
v = 0
cost(1) = 0.0_nag_wp

!         Loop over all possible change point locations
Do i = floc, lloc
this_cost = ruser(2*i-1) + ruser(2*i)

If (this_cost<cost(1) .Or. v==0) Then
!             Update the proposed change point location
v = i
cost(1) = this_cost
cost(2) = ruser(2*i-1)
cost(3) = ruser(2*i)
End If
End Do
Else
!         Need to calculate the cost for the full segment
cost(1) = ruser(2*w-1)
!         No need to populate the rest of COST or V
End If

!       Reset the IEEE exception mode back to what it was
Call x07cbf(cexmode)

!       Set info nonzero to terminate execution for any reason
info = 0
End Subroutine chgpfn

Subroutine get_data(nin,n,y,iuser,ruser)
!       Read in data that is specific to the cost function

!       .. Scalar Arguments ..
Integer, Intent (In)           :: n, nin
!       .. Array Arguments ..
Real (Kind=nag_wp), Allocatable, Intent (Out) :: ruser(:), y(:)
Integer, Allocatable, Intent (Out) :: iuser(:)
!       .. Local Scalars ..
Real (Kind=nag_wp)             :: shape
!       .. Executable Statements ..

!       Read in the series of interest
Allocate (y(1:n))

!       Read in the shape parameter for the Gamma distribution

!       We are going to use RUSER for two purposes, firstly to store the shape
!       parameter, and we also need an additional 2*N elements of workspace
!       we reference from 0 to make the coding easier later

!       IUSER is not going to be used
Allocate (iuser(0),ruser(0:2*n))

!       Store the shape parameter in the 0th element of RUSER
ruser(0) = shape

!       We will be populating the other elements of RUSER in the first
!       call to CHGPFN

Return
End Subroutine get_data
End Module g13nefe_mod

Program g13nefe

!     .. Use Statements ..
Use g13nefe_mod, Only: chgpfn, get_data
Use nag_library, Only: g13nef, nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Parameters ..
Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: beta
Integer                          :: i, ifail, mdepth, minss, n, ntau
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: ruser(:), y(:)
Integer, Allocatable             :: iuser(:), tau(:)
!     .. Intrinsic Procedures ..
Intrinsic                        :: repeat
!     .. Executable Statements ..
Write (nout,*) 'G13NEF Example Program Results'
Write (nout,*)

!     Skip heading in data file

!     Read in the problem size, penalty, minimum segment size and
!     maximum depth
Read (nin,*) n, beta, minss, mdepth

!     Read in the rest of the data, that (may be) dependent on the cost
!     function
Call get_data(nin,n,y,iuser,ruser)

!     Allocate output arrays
Allocate (tau(n))

!     Call routine to detect change points
ifail = 0
Call g13nef(n,beta,minss,mdepth,chgpfn,ntau,tau,y,iuser,ruser,ifail)

!     Display the results
Write (nout,99999) ' -- Change Points --'
Write (nout,99999) ' Number     Position'
Write (nout,99999) repeat('=',21)
Do i = 1, ntau
Write (nout,99998) i, tau(i)
End Do

99999 Format (1X,A)
99998 Format (1X,I4,7X,I6)
End Program g13nefe
```