forums.silverfrost.com

[Mendo910603]

Hi. I'm having troubles with my code. I'm using Numerical Recipes subroutines in my program but at the time of execute the program, I have the following error:

SALFORD RUN-TIME LIBRARY
Stack overflow: Re-link with bigger stack value
Stack:reserve.commit
.will attempt to trace back

I haven't find the error in the code, I don't know if you can help me. Thanks
Here's the code:

PROGRAM Bulirsch_Stoer
Use NRTYPE
INTERFACE
subroutine bsstep(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs)
use nrtype
real(sp),dimension(,intent(inout)::y
real(sp),dimension(,intent(in)::dydx,yscal
real(sp),intent(inout):
real(sp),intent(in)::htry,eps
real(sp),intent(out)::hdid,hnext
end subroutine bsstep
END INTERFACE
PARAMETER(N=4)
real,dimension(n)::y,dydx,yscal
X=1.0
Y(1)=bessj0_s(x)
Y(2)=bessj1_s(x)
Y(3)=bessj_s(2,x)
Y(4)=bessj_s(3,x)
DYDX(1)=-Y(2)
DYDX(2)=Y(1)-Y(2)
DYDX(3)=Y(2)-2.0*Y(3)
DYDX(4)=Y(3)-3.0*Y(4)
DO I=1,4
YSCAL(I)=1.0
CONTINUE
HTRY=1.0
WRITE(*,'(/1X,T8,A,T19,A,T31,A,T43,A/)')
END DO
DO I=1,15
EPS=EXP(-FLOAT(I))
CALL bsstep(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs)
WRITE(*,'(2X,E12.4,F8.2,2X,2F12.6)') EPS,HTRY,HDID,HNEXT
CONTINUE
ENDDO
END PROGRAM Bulirsch_Stoer
SUBROUTINE DERIVS(X,Y,DYDX)
PARAMETER(N=4)
real,dimension(n)::y,dydx
DYDX(1)=-Y(2)
DYDX(2)=Y(1)-(1.0/X)*Y(2)
DYDX(3)=Y(2)-(2.0/X)*Y(3)
DYDX(4)=Y(3)-(3.0/X)*Y(4)
End subroutine derivs
subroutine BSSTEP(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs)
use NRTYPE;use nrutil,only:arth,assert_eq,cumsum,iminloc,nrerror,&
outerdiff,outerprod,upper_triangle
use NR,only:mmid,pzextr
implicit none
real(sp),dimension(,intent(inout)::y
real(sp),dimension(,intent(in)::dydx,yscal
real(sp),intent(inout):
real(sp),intent(in)::htry,eps
real(sp),intent(out)::hdid,hnext
interface
SUBROUTINE derivs(x,y,dydx)
use NRTYPE
implicit none
real(sp),intent(in):
real(sp),dimension(,intent(in)::y
real(sp),dimension(,intent(out)::dydx
end subroutine derivs
end interface
integer(I4B),parameter::imax=9,kmaxx=imax-1
real(sp),parameter::safe1=0.25_sp,safe2=0.7_sp,redmax=1.0e-5_sp,&
redmin=0.7_sp,tiny=1.0e-30_sp,scalmx=0.1_sp
integer(I4B)::k,km,iq,ndum
integer(I4B),dimension(imax)::nseq=(/2,4,6,8,10,12,14,16,18/)
integer(I4B),save::kopt,kmax
real(sp),dimension(kmaxx,kmaxx),save::alf
real(sp),dimension(kmaxx)::err
real(sp),dimension(imax),save::a
real(sp),save::epsold=-1.0_sp,xnew
real(sp)::eps1,errmax,fact,h,red,scale,wrkmin,xest
real(sp),dimension(size(y))::yerr,ysav,yseq
logical(LGT)::reduct
logical(LGT),save::first=.true.
ndum = assert_eq(size(y),size(dydx),size(dydx),size(yscal),'bsstep')
if(eps/=epsold) then
hnext=-1.0e29_sp
xnew=-1.0e29_sp
eps1=safe1*eps
a(=cumsum(nseq,1)
do iq=2,KMAXX
a(k+1)=a(k)+nseq(k+1)
end do
do iq=2,KMAXX
do k=1,iq-1
alf(k,iq)=eps1**((a(k+1)-a(iq-1))/&
((a(iq+1)-a(1)+1.0_sp)*(2*k+1)))
enddo
enddo
epsold=eps
do kopt=2,kmaxx-1
if(a(kopt+1)>a(kopt)*alf(kopt-1,kopt)) exit
end do
kmax=kopt
end if
h=htry
ysav(=y(
if(h/=hnext .or. x/=xnew) then
first=.true.
kopt=kmax
end if
reduct=.false.
main_loop: do
do k=1,kmax
xnew=x+h
if(xnew==x)call nrerror('step size underflow in bsstep')
call mmid(ysav,dydx,x,h,nseq(k),yseq,derivs)
xest=(h/nseq(k))**2
call pzextr(k,xest,yseq,y,yerr)
if(k/=1)then
errmax=maxval(abs(yerr(/yscal())
errmax=max(tiny,errmax)/eps
km=k-1
err(km)=(errmax/safe1)**(1.0_sp/(2*km+1))
end if
if(k/=1.and.(k>=kopt-1.or.first))then
if(errmax<1.0)exit main_loop
if(k == kmax .or. k == kopt+1) then
red=safe2/err(km)
exit
else if(k==kopt)then
if(alf(kopt-1,kopt)<err(km))then
red=1.0_sp/err(km)
exit
end if
else if(kopt==kmax)then
if(alf(km,kmax-1)<err(km))then
red=alf(km,kmax-1)*saf

[JohnCampbell]

The most often cause of stack overflow is large local arrays.
This can be overcome by transferring these arrays from the stack to other memory areas.
Options include:
a) declare arrays in common or in a module, or
b) make large temporary arrays as ALLOCATABLE and include the extra line of code to ALLOCATE the arrays.

If your program is recursive, then you need to check how to manage these arrays.

I have always found the option of changing the stack size as only a temporary fix that soon fails again.

John

[JohnCampbell]

I had a look at the code you have posted and wonder how have you declared derivs as being an external function in PROGRAM Bulirsch_Stoer .

You appear to have declared it as such in subroutine bsstep, but not in the Program or for the bsstep interface ?

I am not familiar with your use of an interface for derivs in BSSTEP and I'm not sure that you need the derivs interface, but suggest you do need an external declaration.

John

[Mendo910603]

The most often cause of stack overflow is large local arrays.
This can be overcome by transferring these arrays from the stack to other memory areas.
Options include:
a) declare arrays in common or in a module, or
b) make large temporary arrays as ALLOCATABLE and include the extra line of code to ALLOCATE the arrays.

If your program is recursive, then you need to check how to manage these arrays.

I have always found the option of changing the stack size as only a temporary fix that soon fails again.

John

I checked my code again and it appears that I have a problem with the array dimensions for y, dydx and yscal in the main program. But now I have problems with a function declared in the BSSTEP subroutine.
Thanks John!

[Mendo910603]

I had a look at the code you have posted and wonder how have you declared derivs as being an external function in PROGRAM Bulirsch_Stoer .

You appear to have declared it as such in subroutine bsstep, but not in the Program or for the bsstep interface ?

I am not familiar with your use of an interface for derivs in BSSTEP and I'm not sure that you need the derivs interface, but suggest you do need an external declaration.

John

I'm not sure if I'm following you, but what you're saying is that I have to declare the derivs subroutine as external or is not necessary because it has been declared in BSSTEP subroutine?

[davidb]

You DON'T need to use External for DERIVS in the main program.

All you need to do is use either External OR an Interface block in the BSTEP subroutine, which you are doing correctly!

However, I couldn't get your code to compile as it is incomplete. If you can post it all, I could take a look. I am reasonably familiar with the NR bstep procedure.
_________________
Programmer in: Fortran 77/95/2003/2008, C, C++ (& OpenMP), java, Python, Perl

[davidb]

One thing I can see is your definition of DERIVS in the Interface block is different from its actual definition. The actual subroutine uses an explicit dimension for Y and DYDX "dimension(n)" but you have used a deferred dimension "dimension(:)" in the interface block. This isn't allowed.

You should be getting a compiler error.
_________________
Programmer in: Fortran 77/95/2003/2008, C, C++ (& OpenMP), java, Python, Perl

[JohnCampbell]

David,

Based on the code sample that has been provided: In program Bulirsch_Stoer, there is an interface definition for bsstep, but this does not contain a definition of derivs.
When bsstep is called (in the program) there is no indication that derivs is an external function.

Subroutine bsstep does declare derivs as a function by way of an interface. I presume that the smile indicates that y and dydx are dimension(:).
As far as I can understand from the code supplied, this information is not available to program Bulirsch_Stoer and would cause the stack error.
It would be a much simpler definition, if derivs used F77 dimension convention, then derivs would only require to be declared external in both routines.
I think a work around would be to include the interface definition for derivs in the main program.

We don't know what is declared in module NRTYPE. It may contain derivs ? Oops, just saw derivs routine so I don't know why an interface is used at all.
You should use external, which could be included in the BSSTEP interface.

John

[davidb]

You don't need to declare DERIVS external in the main program. This is a widely held belief amongst Fortran programmers, especially ones who started with Fortran 77, but its wrong. (Edit: correction. If you don't use external, you need an interface block for DERIVS or it needs to have been declared in a module).

(If derivs was a function you would have to declare its return type).

What is needed is either declare DERIVS External in subroutine BSTEP or provide a valid interface block for it in BSTEP.

Both methods tell the compiler that dummy argument DERIVS is a subroutine. Using External is the simplest way and what most people seem to use (and what the books teach), but an interface block is also allowed instead.

In this case, the interface block supplied inside BSTEP doesn't match the interface in the DERIVS subroutine. The code should not even compile. (Edit: but perhaps only /CHECKMATE will catch this).

When the call is made from BSTEP to DERIVS it must pass an array descriptor (also called a dope vector) indicating the size and shape of Y and DYDX. These "hidden" arguments are pushed on the stack. However, the actual DERIVS subroutine doesn't expect these so cannot manage the stack correctly, probably leading to the stack error.

If only the full code was posted, it could simply be corrected.

BTW. I know what's in NRTYPE (have the book )
_________________
Programmer in: Fortran 77/95/2003/2008, C, C++ (& OpenMP), java, Python, Perl

[JohnCampbell]

David,

When derivs is called from within bsstep, it's external attribute is defined (incorrectly) by the interface, although an EXTERNAL would suffice.

However, when BSSTEP is called from the main program, it provides an argument derivs. If there is no EXTERNAL attribute for derivs defined in the main program, what would the compiler provide ? I would expect it is a real variable, as there is no altenative definition of derivs available.
The INTERFACE provided for BSSTEP in the main program does not include any implied function definition for derivs. This interface is the limit of information provided to the compiler about derivs in the main program and it is not sufficient.

John

[davidb]

Yes, you're right . I missed that.

There needs to be something which says that DERIVS is a subroutine in the main program.

Your approach with external would suffice.

If the OP wants to keep the interface block for DERIVES, then would have to add "External DERIVS" both to the main program and to the interface block for BSTEP.

If the OP wants to use explicit interfaces instead of External, the solution would be to nest an interface block for DERIVS inside the one for BSTEP.

Here is a simple example using nested interface blocks.

[Mendo910603]

Thanks John and Davidb for your comments. I followed what you've told me, and I could solve the problem, but it appeared another one related to a function declared in the nrutil file named assert_eq. This function compares the dimensions of the arrays y,dydx and yscal. If they are the same, then the code follows its path, otherwise the subroutine bsstep calls a function called nrerror and a error message appears.
I've checked the code many times, but I cannot see the error, because either the subroutine or the main program the dimensions of the arrays are the same. I don't know if you can help me because I cannot figure it out.
The error is the following:
nrerror: an assert_eq failed with this tag:bsstep
STOP program terminated by assert_eq3
And the code is:
- Main program
Program D15R6
External Derivs
Parameter(n=4)
dimension y(n),dydx(n),yscal(n)
x=1.0
y(1)=bessj0_s(x)
y(2)=bessj1_s(x)
y(3)=bessj_s(2,x)
y(4)=bessj_s(3,x)
dydx(1)=-y(2)
dydx(2)=y(1)-y(2)
dydx(3)=y(2)-2.0*y(3)
dydx(4)=y(3)-3.0*y(4)
do i=1,n
yscal(i)=1.0
continue
end do
htry=1.0
write(*,'(/1X,T8,A,T19,A,T31,A,T43,A/)')
Do i=1,15
eps=exp(float(i))
call bsstep(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs)
write(*,'(2X,E12.4,F8.2,2X,2F12.6)') EPS,HTRY,HDID,HNEXT
Continue
end do
End program D15R6
Subroutine derivs(x,y,dydx)
dimension y(4),dydx(4)
dydx(1)=-y(2)
dydx(2)=y(1)-(1.0/x)*y(2)
dydx(3)=y(2)-(2.0/x)*y(3)
dydx(4)=y(3)-(3.0/x)*y(4)
return
end subroutine derivs

-Subroutine BSSTEP
Subroutine bsstep(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs)
Use nrtype;Use nrutil,ONLY:arth,assert_eq,cumsum,iminloc,nrerror, &
outerdiff,outerprod,upper_triangle
Use nr,only:mmid,pzextr
Implicit none
Real(sp),dimension(,intent(inout)::y
Real(sp),dimension(,intent(in)::dydx,yscal
Real(sp),intent(inout):
Real(sp),intent(in)::htry,eps
Real(sp),intent(out)::hdid,hnext
Interface
Subroutine derivs(x,y,dydx)
Use nrtype
Implicit none
Real(sp),intent(in):
Real(sp),dimension(,intent(in)::y
Real(sp),dimension(,intent(out)::dydx
End Subroutine derivs
End Interface
Integer(I4B),parameter::IMAX=9,KMAXX=IMAX-1
Real(sp),parameter::SAFE1=0.25_sp,SAFE2=0.7_sp,REDMAX=1.0e-5_sp,&
REDMIN=0.7_sp,TINY=1.0e-30_sp,SCALMX=0.1_sp
Integer(I4B)::k,km,ndum
Integer(I4B),dimension(IMAX)::nseq=(/ 2,4,6,8,10,12,14,16,18 /)
Integer(I4B),save::kopt,kmax
Real(sp),dimension(KMAXX,KMAXX),save::alf
Real(sp),dimension(KMAXX)::err
Real(sp),dimension(IMAX),save::a
Real(sp),save::epsold=1.0_sp,xnew
Real(sp)::eps1,errmax,fact,h,red,scale,wrkmin,xest
Real(sp),dimension(size(y))::yerr,ysav,yseq
Logical(LGT)::reduct
Logical(LGT),save::first=.true.
ndum=assert_eq(size(y),size(dydx),size(yscal),'bsstep')
if (eps/=epsold) then
hnext=-1.0e-29_sp
xnew=-1.0e-29_sp
eps1=SAFE1*eps
a(=cumsum(nseq,1)
where (upper_triangle(KMAXX,KMAXX)) alf=eps1** &
(outerdiff(a(2,a(2)/outerprod(arth( &
3.0_sp,2.0_sp,KMAXX),(a(2-a(1)+1.0_sp)))
epsold=eps
do kopt=2,KMAXX-1
if (a(kopt+1)>a(kopt)*alf(kopt-1,kopt)) exit
end do
kmax=kopt
end if
h=htry
ysav(=y(
if (h/=hnext.or.x/=xnew) then
first=.true.
kopt=kmax
end if
reduct=.false.
main_loop: do
do k=1,kmax
xnew=x+h
if (xnew==x) call nrerror('step size underflow in bsstep')
call mmid(ysav,dydx,x,h,nseq(k),yseq,derivs)
xest=(h/nseq(k))**2
call pzextr(k,xest,yseq,y,yerr)
if (k/=1) then
errmax=maxval(abs(yerr(/yscal())
errmax=max(TINY,errmax)/eps

[Mendo910603]

km=k-1
err(km)=(errmax/SAFE1)**(1.0_sp/(2*km+1))
end if
if (k/=1 .and. (k>=kopt-1 .or. first)) then
if (errmax < 1.0) exit main_loop
if (k==kmax .or. k==kopt+1) then
red=SAFE2/err(km)
exit
else if (k==kopt) then
if (alf(kopt-1,kopt) < err(km)) then
red=1.0_sp/err(km)
exit
end if
else if (kopt==kmax) then
if (alf(km,kmax-1) < err(km)) then
red=alf(km,kmax-1)*SAFE2/err(km)
exit
end if
else if (alf(km,kopt) < err(km)) then
red=alf(km,kopt-1)/err(km)
exit
end if
end if
end do
red=max(min(red,REDMIN),REDMAX)
h=h*red
reduct=.true.
end do main_loop
x=xnew
hdid=h
first=.false.
kopt=1+iminloc(a(2:km+1)*max(err(1:km),SCALMX))
scale=max(err(kopt-1),SCALMX)
wrkmin=scale*a(kopt)
hnext=h/scale
if (kopt >= k .and. kopt /= kmax .and. .not. reduct) then
fact=max(scale/alf(kopt-1,kopt),SCALMX)
if (a(kopt+1)*fact <= wrkmin) then
hnext=h/fact
kopt=kopt+1
end if
end if
End subroutine bsstep

-Function assert_eq contained in the nrutil file and the one who launches the error message
FUNCTION assert_eq3(n1,n2,n3,string)
CHARACTER(LEN=*), INTENT(IN) :: string
INTEGER, INTENT(IN) :: n1,n2,n3
INTEGER :: assert_eq3
if (n1 == n2 .and. n2 == n3) then
assert_eq3=n1
else
write (*,*) 'nrerror: an assert_eq failed with this tag:', &
string
STOP 'program terminated by assert_eq3'
end if
END FUNCTION assert_eq3

[Mendo910603]

Thanks for your comments and your help

[JohnCampbell]

A problem is that the array sizes are not available for the call to assert_eq3.
I think the problem is that in BSSTEP you have retained the definition of y, dxdy and yscal as dimension(:), but not provided an interface for the BSSTEP call in the main program.
My recommendation is to bring back the INTERFACE for BSSTEP, but convert the INTERFACE for DERIVS to an EXTERNAL, as dimension(4) (in derivs) does not require an interface.
Assuming the USEd modules do not relate to this, the code might look like: