Topic: QUADPACK in SF => SIGSEGV error

Has anyone tried using QUADPACK integration routines in SF?

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

  • jeff
  • Administrator
  • Offline

Re: QUADPACK in SF => SIGSEGV error

Looking at the definition of dqag, it looks like work should be an array rather than a scalar as you've defined it. The invalid memory reference is probably a result of this.   Additionally, iwork should be an integer array as well.

Jeff Armstrong
Approximatrix, LLC

Re: QUADPACK in SF => SIGSEGV error

Thank you for pointing out the stupid mistake by me. Correcting it, however, did not cure the segmentation fault problem. The modified main program is below.

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

  • jeff
  • Administrator
  • Offline

Re: QUADPACK in SF => SIGSEGV error

It looks like your array dimensions are switched.  You should have:

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.

Jeff Armstrong
Approximatrix, LLC

Re: QUADPACK in SF => SIGSEGV error

Thank you very much for your help. My experience is mainly from Matlab and Fortran is far more complicated, especially debugging.

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...

  • jeff
  • Administrator
  • Offline

Re: QUADPACK in SF => SIGSEGV error

There isn't really a mechanism for anonymous functions ala MATLAB in Fortran.  You'd have to use public variables.  My usual method for getting around this limitation is to use a module that contains data, set the values in this module prior to calling the integrator, then finally "use"-ing the module in the function to be integrated, providing access to the data.  This method keeps the data contained strictly in the module.  It might not be ideal, but it might help.

Jeff Armstrong
Approximatrix, LLC