I have one quite general problem concerning external routines, such as integration routines. For example, I want to integrate function f(x,a,b,c,...) over x. Routines like dqag want a function f(x). I have not yet found any good way to tell the function f(x) values for the parameters a,b,c,... Public variables are naturally one option, but are they really the best solution? I should calculate lots of integrals with different values of a,b,c...

lenw = 400
limit = 100
allocate(iwork(limit))
allocate(work(lenw))

The SIGSEGV is being called because the integration routines were accessing values outside of work's bounds.  It only survives by chance.

Sometimes the program manages to calculate the integral, and then the segmentation fault occurs AFTER I have answered the question "continue program?". If I comment the "call dqag" and "write" lines, then there are no errors. I am a beginner in Fortran, and I do not understand how can it be possible that a program fails after it has executed all the lines? The interface section does not seem to make any difference.

The dqag instructions say that "testfunction" should be declared as "external", but that does not seem to be possible, because that is a module function.

program quadpacktest
use testmodule
use constants
implicit none
!interface
!    subroutine dqag(testfunction,a,b,epsabs,epsrel,key,y,abserr,neval,ier,limit,lenw,last,iwork,work)
!        use constants
!        implicit none
!        real (dp), external :: testfunction
!        real (dp), intent(in) :: a,b,epsabs,epsrel
!        integer, intent(in) :: key, limit, lenw
!        real (dp), intent(out) :: y, abserr
!        integer, intent(out) :: neval, ier, last
!        real (dp), intent(out), allocatable :: work(:)
!        integer, intent(out), allocatable :: iwork(:)
!    end subroutine dqag
!end interface
real (dp) :: x, y
real (dp) :: a, b, epsabs, epsrel, abserr
real (dp), allocatable :: work(:)
integer :: key, neval, ier, limit, lenw, last
integer, allocatable :: iwork(:)
character :: line
x=2.0_dp
a=0.0_dp
b=1.0_dp
epsabs=1.0e-3_dp
epsrel=1.0e-3_dp
key=2
limit=100
allocate(work(limit))
lenw=400
allocate(iwork(lenw))
!These two lines I sometimes comment:
call dqag(testfunction,a,b,epsabs,epsrel,key,y,abserr,neval,ier,limit,lenw,last,iwork,work)
write(*,*) y, abserr, neval, ier
DO
WRITE (*,*) 'Continue program (Y,y)?'
READ (*,'(A)') line
line = ADJUSTL(line)
IF (line(1:1) == 'Y' .OR. line(1:1) == 'y') EXIT
END DO
end program quadpacktest

I wrote a simple test program trying to make QUADPACK work. After locating plenty of additional subroutines (for example, reference BLAS and XERROR from Netlib), my program finally compiles okay, but when I try to run it, I obtain the following error message:

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
#0  ffffffff

My program and modules are below:

program quadpacktest
use testmodule
use constants
implicit none
real (dp) :: x, y
real (dp) :: a, b, epsabs, epsrel, abserr, ier, work
integer :: neval, key, limit, lenw, last, iwork
x=2.0_dp
a=0.0_dp
b=1.0_dp
epsabs=1.0e-3_dp
epsrel=1.0e-3_dp
key=2
limit=100
lenw=400
iwork=100
call dqag(testfunction,a,b,epsabs,epsrel,key,y,abserr,neval,ier,limit,lenw,last,iwork,work)
write(*,*) y
end program quadpacktest

module testmodule
use constants
implicit none
contains
pure function testfunction(x) result(y)
    real (dp), intent(in) :: x
    real (dp) :: y
    y=x**1.5_dp
end function testfunction
end module testmodule

module constants
implicit none
integer, parameter, public :: dp = kind(0.d0)
real (dp), parameter, public :: pi = 3.1415926535897932_dp, e=1.602176462e-19_dp, kb=1.3806503e-23_dp
end module constants