NAG Library Manual, Mark 28.3
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
!   Mark 28.3 Release. NAG Copyright 2022.

    Module e04gnfe_mod
!     .. Use Statements ..
      Use iso_c_binding, Only: c_ptr
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: lsqfun, lsqgrd

    Contains

      Subroutine lsqfun(nvar,x,nres,rx,inform,iuser,ruser,cpuser)
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (In)      :: cpuser
        Integer, Intent (Inout)        :: inform
        Integer, Intent (In)           :: nres, nvar
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (Out) :: rx(nres)
        Real (Kind=nag_wp), Intent (In) :: x(nvar)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Integer                        :: i
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sin
!       .. Executable Statements ..

        rx(:) = 0.0_nag_wp

        Do i = 1, nres
          rx(i) = ruser(nres+i) - x(1)*ruser(i)*sin(-x(2)*ruser(i))
        End Do

        inform = 0
      End Subroutine lsqfun

      Subroutine lsqgrd(nvar,x,nres,nnzrd,rdx,inform,iuser,ruser,cpuser)
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (In)      :: cpuser
        Integer, Intent (Inout)        :: inform
        Integer, Intent (In)           :: nnzrd, nres, nvar
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: rdx(nnzrd), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: x(nvar)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Integer                        :: i
!       .. Intrinsic Procedures ..
        Intrinsic                      :: cos, sin
!       .. Executable Statements ..

        rdx(:) = 0.0_nag_wp

        Do i = 1, nres
          rdx((i-1)*nvar+1) = -ruser(i)*sin(-x(2)*ruser(i))
          rdx((i-1)*nvar+2) = ruser(i)**2*x(1)*cos(x(2)*ruser(i))
        End Do

        inform = 0
      End Subroutine lsqgrd

    End Module e04gnfe_mod

    Program e04gnfe

!     .. Use Statements ..
      Use e04gnfe_mod, Only: lsqfun, lsqgrd
      Use iso_c_binding, Only: c_null_ptr, c_ptr
      Use nag_library, Only: e04gnf, e04gnu, e04gnx, e04gny, e04raf, e04rhf,   &
                             e04rjf, e04rmf, e04rzf, e04zmf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: infbnd = 1.0E20_nag_wp
      Integer, Parameter               :: nout = 6
!     .. Local Scalars ..
      Type (c_ptr)                     :: cpuser, handle
      Integer                          :: idlc, ifail, isparse, nclin, nnza,   &
                                          nnzrd, nres, nvar, x_idx
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:), bla(:), blx(:), bua(:),        &
                                          bux(:), ruser(:), rx(:), x(:)
      Real (Kind=nag_wp)               :: rinfo(100), stats(100)
      Integer, Allocatable             :: icola(:), irowa(:)
      Integer                          :: icolrd(0), irowrd(0), iuser(0)
!     .. Executable Statements ..

      Write (nout,*) 'e04gnf Example Program Results'
      Write (nout,*)
      Flush (nout)

      cpuser = c_null_ptr

!     Problem size
      nvar = 2
!     Residual quantity
      nres = 24
      Allocate (ruser(2*nres))
!     Data for model y = f(t;A,w) := A*t*sin(-w*t)
!     Data generated using A = 0.1 and w = 0.8 and added random noise
!     Residual 4, 12, 16 and 20 are outliers
!     t(:) =
      ruser(1:nres) = (/1.0E+00_nag_wp,2.0E+00_nag_wp,3.0E+00_nag_wp,          &
        4.0E+00_nag_wp,5.0E+00_nag_wp,6.0E+00_nag_wp,7.0E+00_nag_wp,           &
        8.0E+00_nag_wp,9.0E+00_nag_wp,1.0E+01_nag_wp,1.1E+01_nag_wp,           &
        1.2E+01_nag_wp,1.3E+01_nag_wp,1.4E+01_nag_wp,1.5E+01_nag_wp,           &
        1.6E+01_nag_wp,1.7E+01_nag_wp,1.8E+01_nag_wp,1.9E+01_nag_wp,           &
        2.0E+01_nag_wp,2.1E+01_nag_wp,2.2E+01_nag_wp,2.3E+01_nag_wp,           &
        2.4E+01_nag_wp/)
!     y(:) =
      ruser(nres+1:nres*2) = (/0.0523_nag_wp,-0.1442_nag_wp,-0.0422_nag_wp,    &
        1.8106_nag_wp,0.3271_nag_wp,0.4684_nag_wp,0.4593_nag_wp,               &
        -0.0169_nag_wp,-0.7811_nag_wp,-1.1356_nag_wp,-0.5343_nag_wp,           &
        -3.0043_nag_wp,1.1832_nag_wp,1.5153_nag_wp,0.7120_nag_wp,              &
        -2.2923_nag_wp,-1.4871_nag_wp,-1.7083_nag_wp,-0.9936_nag_wp,           &
        -5.2873_nag_wp,1.7555_nag_wp,2.0642_nag_wp,0.9499_nag_wp,              &
        -0.6234_nag_wp/)

      iuser(:) = 0

!     Initialize handle
      ifail = 0
      Call e04raf(handle,nvar,ifail)

!     Define residuals structure, isparse=0 means the residual structure is
!     dense => irowrd and icolrd are not accessed
      isparse = 0
      nnzrd = 0
      ifail = 0
      Call e04rmf(handle,nres,isparse,nnzrd,irowrd,icolrd,ifail)

!     Set options: use least-squares loss function
      Call e04zmf(handle,'NLDF Loss Function Type = L2',ifail)
!     Change printed output verbosity
      Call e04zmf(handle,'Print Level = 1',ifail)
      Call e04zmf(handle,'Print Options = No',ifail)

!     Define starting point
      Allocate (x(nvar),rx(nres))
      x(1:nvar) = (/0.3_nag_wp,0.7_nag_wp/)

!     Define bounds
      Allocate (blx(nvar),bux(nvar))
      blx(1) = -1.0_nag_wp
      bux(1) = infbnd
      blx(2) = 0.0_nag_wp
      bux(2) = 1.0_nag_wp
      ifail = 0
      Call e04rhf(handle,nvar,blx,bux,ifail)

!     Add linear constraints
      nclin = 1
      nnza = 2
      Allocate (bla(nclin),bua(nclin),irowa(nnza),icola(nnza),a(nnza))
      bla(1) = 0.2_nag_wp
      bua(1) = 1.0_nag_wp
      irowa(1:nnza) = (/1,1/)
      icola(1:nnza) = (/1,2/)
      a(1:nnza) = (/1.0_nag_wp,1.0_nag_wp/)
      idlc = 0
      ifail = 0
      Call e04rjf(handle,nclin,bla,bua,nnza,irowa,icola,a,idlc,ifail)

      Write (nout,*)                                                           &
        'First solve the problem using least squares loss function'
      Write (nout,*)                                                           &
        '---------------------------------------------------------'
!     Call the solver
      ifail = -1
      Call e04gnf(handle,lsqfun,lsqgrd,e04gnx,e04gny,e04gnu,nvar,x,nres,rx,    &
        rinfo,stats,iuser,ruser,cpuser,ifail)

!     Print solution if optimal or suboptimal solution found
      If (ifail==0 .Or. ifail==50) Then
        Write (nout,*)
        Write (nout,99999) 'Optimal parameters X (Least Squares):'
        Write (nout,99997) 'x_idx', ' Value '
        Do x_idx = 1, nvar
          Write (nout,99998) x_idx, x(x_idx)
        End Do
      End If

      Write (nout,*)
      Write (nout,*) 'Now solve the problem using SmoothL1 loss function'
      Write (nout,*)                                                           &
        '---------------------------------------------------------'
!     Use robust nonlinear regression and resolve the model
      Call e04zmf(handle,'NLDF Loss Function Type = SmoothL1',ifail)

!     Call the solver to resolve
      x(1:nvar) = (/0.3_nag_wp,0.7_nag_wp/)
      ifail = -1
      Call e04gnf(handle,lsqfun,lsqgrd,e04gnx,e04gny,e04gnu,nvar,x,nres,rx,    &
        rinfo,stats,iuser,ruser,cpuser,ifail)

!     Print solution if optimal or suboptimal solution found
      If (ifail==0 .Or. ifail==50) Then
        Write (nout,*)
        Write (nout,99999) 'Optimal parameters X (SmoothL1):'
        Write (nout,99997) 'x_idx', ' Value '
        Do x_idx = 1, nvar
          Write (nout,99998) x_idx, x(x_idx)
        End Do
      End If

!     Free the handle memory
      ifail = 0
      Call e04rzf(handle,ifail)

99999 Format (1X,A)
99998 Format (2X,I5,3X,Es7.2e1)
99997 Format (2X,A5,3X,A7)

    End Program e04gnfe