Stack overflow problem

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)::x 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.0Y(3) DYDX(4)=Y(3)-3.0Y(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)::x 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)x real(sp),dimension(:),intent(in)y real(sp),dimension(:),intent(out)dydx end subroutine derivs end interface integer(I4B),parameterimax=9,kmaxx=imax-1 real(sp),parametersafe1=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),savekopt,kmax real(sp),dimension(kmaxx,kmaxx),savealf real(sp),dimension(kmaxx)err real(sp),dimension(imax),savea real(sp),saveepsold=-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),savefirst=.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=safe1eps 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)(2k+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/(2km+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(kkopt)then if(alf(kopt-1,kopt)<err(km))then red=1.0_sp/err(km) exit end if else if(koptkmax)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=hred 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=scalea(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 SUBROUTINE mmid(y,dydx,xs,htot,nstep,yout,derivs) USE nrtype; USE nrutil, ONLY: assert_eq3,swap IMPLICIT NONE INTEGER(I4B),INTENT(IN)::nstep REAL(SP),INTENT(IN)::xs,htot REAL(SP),DIMENSION(:),INTENT(IN)::y,dydx REAL(SP),DIMENSION(:),INTENT(OUT)::yout INTERFACE SUBROUTINE derivs(x,y,dydx) USE nrtype IMPLICIT NONE REAL(SP),INTENT(IN)::x REAL(SP),DIMENSION(:),INTENT(IN)::y REAL(SP),DIMENSION(:),INTENT(OUT)::dydx END SUBROUTINE derivs END INTERFACE INTEGER(I4B)::n,ndum REAL(SP)::h,h2,x REAL(SP),DIMENSION(size(y))::ym,yn ndum=assert_eq3(size(y),size(dydx),size(yout),'mmid') h=htot/nstep ym=y yn=y+hdydx x=xs+h call derivs(x,yn,yout) h2=2.0_sph do n=2,nstep call swap(ym,yn) yn=yn+h2yout x=x+h call derivs(x,yn,yout) enddo yout=0.5_sp*(ym+yn+hyout) END SUBROUTINE mmid SUBROUTINE pzextr(iest,xest,yest,yz,dy) USE nrtype; USE nrutil, ONLY : assert_eq,nrerror IMPLICIT NONE INTEGER(I4B),INTENT(IN)::iest REAL(SP),INTENT(IN)::xest REAL(SP),DIMENSION(:),INTENT(IN)::yest REAL(SP),DIMENSION(:),INTENT(OUT)yz,dy INTEGER(I4B),PARAMETERIEST_MAX=16 INTEGER(I4B)j,nv INTEGER(I4B),SAVEnvold=-1 REAL(SP)delta,f1,f2 REAL(SP),DIMENSION(size(yz))d,tmp,q REAL(SP),DIMENSION(IEST_MAX),SAVEx REAL(SP),DIMENSION(:,:),ALLOCATABLE,SAVEqcol nv=assert_eq(size(yz),size(yest),size(dy),'pzextr') if (iest>IEST_MAX) call & nrerror('pzextr:probable misuse, too much extrapolation') if (nv/=nvold) then if (allocated(qcol))deallocate(qcol) allocate(qcol(nv,IEST_MAX)) nvold=nv end if x(iest)=xest dy(:)=yest(:) yz(:)=yest(:) if (iest==1) then qcol(:,1)=yest(:) else d(:)=yest(:) do j=1,iest-1 delta=1.0_sp/(x(iest-j)-xest) f1=xestdelta f2=x(iest-j)delta q(:)=qcol(:,j) qcol(:,j)=dy(:) tmp(:)=d(:)-q(:) dy(:)=f1tmp(:) d(:)=f2*tmp(:) yz(:)=yz(:)+dy(:) end do qcol(:,iest)=dy(:) end if END SUBROUTINE pzextr

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:

1. declare arrays in common or in a module, or
2. 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 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

*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:

1. declare arrays in common or in a module, or
2. 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!

*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?

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

All you need to do is use either External [u:c2864f4939]OR[/u:c2864f4939] 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.

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.

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

You [u:2a0471e1ca]don't[/u:2a0471e1ca] 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 [u:2a0471e1ca]is[/u:2a0471e1ca] needed is either declare DERIVS External in subroutine BSTEP or provide a [u:2a0471e1ca]valid[/u:2a0471e1ca] 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 😃 )

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

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.

program anon

   interface
      subroutine my_abs(x, y)
         real x, y
      end subroutine my_abs
   end interface

   interface
      subroutine eval_func(x, y, fnc)
         real, intent(in) :: x
         real, intent(out) :: y
         interface
            subroutine fnc(x, y)
               real x, y
            end subroutine fnc
         end interface
      end subroutine eval_func
   end interface
   
   call eval_func(-1.0, y, my_abs)
   
   print *,'y = ', y, ' (should be 1.0)'
   
end program

! Subroutine to evaluate a function (provided as subroutine)
subroutine eval_func(x, y, fnc)
   real, intent(in) :: x
   real, intent(out) :: y
   interface
      subroutine fnc(x, y)
         real x, y
      end subroutine fnc
   end interface
   
   call fnc(x, y)
end subroutine eval_func

! Example subroutine passed to eval_func
subroutine my_abs(x, y)
   real x, y
   y = abs(x)
end subroutine my_abs

Using external only is obviously less verbose!

Fortran 2003 has added an 'import' statement which avoids the need for the duplication in the above code, but this isn't supported as an extension in FTN95. An alternative approach, which I would prefer, would be to use modules.

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.0y(3) dydx(4)=y(3)-3.0y(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)::x 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)x Real(sp),dimension(:),intent(in)y Real(sp),dimension(:),intent(out)dydx End Subroutine derivs End Interface Integer(I4B),parameterIMAX=9,KMAXX=IMAX-1 Real(sp),parameterSAFE1=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),savekopt,kmax Real(sp),dimension(KMAXX,KMAXX),savealf Real(sp),dimension(KMAXX)err Real(sp),dimension(IMAX),savea Real(sp),saveepsold=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),savefirst=.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=SAFE1eps 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 km=k-1 err(km)=(errmax/SAFE1)(1.0_sp/(2km+1)) end if if (k/=1 .and. (k>=kopt-1 .or. first)) then if (errmax < 1.0) exit main_loop if (kkmax .or. kkopt+1) then red=SAFE2/err(km) exit else if (kkopt) then if (alf(kopt-1,kopt) < err(km)) then red=1.0_sp/err(km) exit end if else if (koptkmax) 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=hred 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=scalea(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 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*

Thanks for all 😃

•        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=scalea(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*

Thanks for your comments and your help 😃

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:

Program D15R6 
Parameter(n=4) 
dimension y(n),dydx(n),yscal(n) 
!
Interface 
 Subroutine bsstep (y,dydx,x,htry,eps,yscal,hdid,hnext,derivs) 
  Real(sp),dimension(:),intent(inout) :: y 
  Real(sp),dimension(:),intent(in)    :: dydx,yscal 
  Real(sp),intent(inout)              :: x 
  Real(sp),intent(in)                 :: htry,eps 
  Real(sp),intent(out)                :: hdid,hnext 
  external derivs
 end Subroutine bsstep 
End Interface 
!
x=1.0 
...
End program D15R6 

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)              :: x 
Real(sp),intent(in)                 :: htry,eps 
Real(sp),intent(out)                :: hdid,hnext 
external derivs
!
Integer(I4B),parameter :: IMAX  = 9
Integer(I4B),parameter :: KMAXX = IMAX-1 

It appears that in BSSTEP, you have coded some generality into the dimensions of y, dxdx and yscal, but then declare them as dimension(4) in derivs. Either this generality is unnecessary or you have to change the code for derivs and also associated INTERFACE. My recomendation would be to adopt F77 syntax for DERIVS and simplify the code as: Subroutine derivs (x,y,dydx,n) integer n dimension y(n),dydx(n) ! 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

This form would only require an external definition and not the interface for the dimensions. (Can 'n' be other than 4 ? if not, why all the extra confusion?)

John

Yes, 'n' can be other than 4, but I'm using that value because I'm introducing 4 differential equations.

Then if n can be different to 4, and in any case; you need to be consistent with the way you handle the declarations for the arrays you are using.

If you declare an array as dimension(:), then you need an interface for this routine, where it is being called, to provide the dimension. If you are consistent in the way you manage dimension(:) it will work ok.

However, if you declare an array as dimension(n), then you should preferably provide the n as an argument, or via other means. This is a simpler approach.

John