An Interface Issue in ELSUNC

I host an archive of Alan Miller’s Fortran Software and recently Professor N. Shamsundar wrote to me to regarding an interface bug in the elsunc.f90 code for for non-linear least squares with upper and lower bounds on the parameters. As I’ve mentioned previously, I have no plans to actively maintain this vast library but I’ll make an effort to fix “obvious” things and this is a pretty clear case.

This excerpt from Prof. Shamsundar’s original email summarizes the issue:

The subroutine for the object function is declared with assumed shape array arguments in d_elsunc.f90. However, the numerous places from which the subroutine is called do not provide an explicit interface, as is required by the Fortran standard for subroutines with assumed shape arguments. Because of this bug, the test problem code builds and runs but aborts with segfaults. The fix is to provide an explicit interface wherever one is called. For example, with the short module code below, one can add USE FUNINT and replace EXTERNAL FFUNC by procedure(optfun) :: ffunc in subroutines LSUNC, NEWUNC, JACDIF, STEPUC, EVREUC and FSUMSQ.

module funint
    abstract interface
       SUBROUTINE optfun(x, n, f, m, ctrl, c, mdc)
          IMPLICIT NONE
          INTEGER, PARAMETER       :: dp = SELECTED_REAL_KIND(12, 60)
          REAL (dp), INTENT(IN)    :: x(:)
          INTEGER, INTENT(IN)      :: n
          REAL (dp), INTENT(OUT)   :: f(:)
          INTEGER, INTENT(IN)      :: m
          INTEGER, INTENT(IN OUT)  :: ctrl
          REAL (dp), INTENT(OUT)   :: c(:,:)
          INTEGER, INTENT(IN)      :: mdc
       END SUBROUTINE optfun
    end interface
end module

This is indeed correct, so I’ll briefly describe the patch I applied. First, note that elsunc.f90 contains the line

USE precision, ONLY : dp

However, no precision module is defined. Therefore, you’ll have to create your own to tell the routines what double precision should mean:

module precision
  implicit none
  integer, parameter :: dp = selected_real_kind(12, 60)
end module precision

But then, if you try to compile and run the example you’ll see an error like the following, due to the lack of an interface for ffunc:

% gfortran -o d_elsunc precision.f90 elsunc.f90 d_elsunc.f90
% ./d_elsunc
How many iterations?
10

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0  0x10b610fbe
#1  0x10b6116d4
#2  0x7fff8ec6e8e9
#3  0x10b6084a1
#4  0x10b5f61c4
#5  0x10b604eb1
#6  0x10b60684d
#7  0x10b60800e
#8  0x10b6082b7
zsh: segmentation fault  ./d_elsunc

The remedy suggested in the email above uses the abstract interface statement from Fortran 2003 and requires adding an additional module. To fix this using only standard Fortran 90 constructs, one can simply replace all occurrences of EXTERNAL ffunc with an explicit interface declaration like this:

interface
  subroutine ffunc(x, n, f, m, ctrl, c, mdc)
    implicit none
    integer, parameter :: dp = selected_real_kind(12, 60)
    real (dp), intent(in) :: x(:)
    integer, intent(in) :: n
    real (dp), intent(out) :: f(:)
    integer, intent(in) :: m
    integer, intent(in out) :: ctrl
    real (dp), intent(out) :: c(:,:)
    integer, intent(in) :: mdc
  end subroutine ffunc
end interface

I have updated elsunc.f90 as described, replacing all instances of EXTERNAL ffunc with the explicit interface above.

It would be cleaner and more general to avoid repeatedly defining the real kind parameter dp and instead add use precision, only: dp to the interface. Any user-provided ffunc functions would also have to be modified to use dp from the precision module in the same way. However, to keep the changes minimal and avoid breaking any existing code (including the test program), I only made the change described above in which dp is repeatedly defined in the interface each time.