Subject:        Re: stacklimits
   Date:        Wed, 23 May 2001 11:08:50 -0400 (EDT)
  From:        Doug Mccune 

Lew,

It is worth making a fortran-90 comment in this context.  This is
because it is my experience that many stack limit problems arise
from "naive" use of fortran-90 features.

Fortran-90 has a very nice tempting feature-- automatic arrays dimensioned
inside a subroutine using calling arguments-- but it is dangerous because
the automatic arrays count against the "stack limit" which is fairly small
on many machines.

A permanent fix, though a bit painful programming-wise, is to turn the
automatic arrays into explicitly allocated/deallocated arrays.  Such
arrays go on the "heap" and you have essentially the whole virtual
address space of the machine at your disposal-- i.e. only swap space
limits the cumulative size.

I learned this the hard way writing and testing fortran-90 3d spline 
routines.  The spline coefficient evaluators looked something like this:

  subroutine mksplin3(f,nx,ny,nz,...)
  real f(8,nx,ny,nz)  ! passed data & compact spline coefficients array
  real wk(64,nx,ny,nz) ! automatic workspace -- local, not passed
    ...

  end subroutine mksplin3

**but** depending on the size of nx,ny,nz, I would easily exceed the stack
limit on various machines.  Hence this was highly non-portable code.  The
fix was simply to do:

  subroutine mksplin3(f,nx,ny,nz,...)
  real f(8,nx,ny,nz)  ! passed
  real, dimension(:,:,:,:), allocatable :: wk

  ...

  allocate(wk(64,nx,ny,nz))   ! allocate when needed

  ...

  deallocate(wk)              ! be sure to deallocate when done

  ...

  end subroutine mksplin3

This is more code, and requires the programmer to explicitly allocate, 
keep track of, and deallocate the memory (failure to deallocate can lead
to serious problems).  But it is safe against stack limits.

fyi ------------ Doug


PS-- in code with complicated error handling logic, it can be complicated
to keep track of whether various arrays have been allocated or not.  One
trick for coping is:

  subroutine mksplin3(...)
  ...

  do
    call foosub(...,ierr)
    if(ierr.ne.0) exit               ! exit do/enddo region on error
    allocate(wk(64,nx,ny,nz))
    ...
    call foosub2(...,ierr)
    if(ierr.ne.0) exit               ! exit do/enddo region on error
    allocate(wk2(nx,ny,nz))
    ...
    exit                             ! exit do/enddo region when done
  enddo

  deallocate(wk,stat=istat)          ! fetch but ignore status code
  deallocate(wk2,stat=istat)         ! fetch but ignore status code

  end subroutine mksplin3

The "stat=" clause renders the trailing deallocate statements harmless if 
the array has not been allocated.  (Many machines will give an error if
you try to deallocate an array that was never allocated in the first
place).  This "stat=" clause method also works for "pointer" 
allocated arrays-- for which the "allocated" logical intrinsic is unsafe
if the pointer was never intitialized.  If there are no pointer-allocated
arrays or if lack of initialization is not a hazard, one can also code 

  if(allocated(wk)) deallocate(wk)

which makes the intent more transparent.

fyi ---------- Doug

PPS-- why compiler writers make life difficult by putting the automatic
arrays on the stack instead of the heap, I don't understand.  But it is
common practice and we have to adjust to it.