Fortran 90中的堆栈溢出

Question

我已经用Fortran 90.一直致力于精美的相当长一段时间写了一个相当大的程序,但今天我上前一步了一个档次,增加问题大小(它是一家集科研非标准FE-求解,如果帮助任何人...)现在我得到"堆栈溢出"错误消息,自然程序终止,而没有给我任何有用的工作.

该程序首先设置所有相关的数组和矩阵,然后完成后,它会将一些关于此的统计信息打印到日志文件中.即使我有一个新的,更大的问题,这个工作正常(尽管有点慢),但随着"数字运算"开始失败.

让我感到困惑的是,那一点上的所有东西都已经分配了(而且没有错误).我不完全确定堆栈是什么(维基百科和这里的几个步骤没有做太多,因为我对计算机的"幕后"工作只有非常基本的知识).

假设我有一些数组初始化为:

INTEGER,DIMENSION(64) :: IA
REAL(8),DIMENSION(:,:),ALLOCATABLE :: AA, BB

在一些初始化例程(即从文件读取输入等)之后被分配为(我存储一些大小整数以便更容易地传递到固定大小的IA中的子例程):

ALLOCATE( AA(N1,N2) , BB(N1,N2) )
IA(1) = N1
IA(2) = N2

这基本上是在初始部分发生的事情,到目前为止一直很好.但是当我接下来调用子程序时

CALL ROUTINE_ONE(AA,BB,IA)

例程看起来像(没什么特别的):

SUBROUTINE ROUTINE_ONE(AA,BB,IA)
IMPLICIT NONE
INTEGER,DIMENSION(64) :: IA
REAL(8),DIMENSION(IA(1),IA(2)) :: AA, BB
...
do lots of other stuff
...
END SUBROUTINE ROUTINE_ONE

现在我收到一个错误!屏幕输出显示:

forrtl: severe (170): Program Exception - stack overflow

然而,当我与调试器中运行该程序它一个名为处断裂线419 winsig.c(不是我的文件,但可能是编译器的一部分？).它似乎是一个sigreterror:被调用的例程的一部分,它是已被调用的默认情况,返回文本Invalid signal or error.这附有评论专栏,奇怪地说/* should never happen, but compiler can't tell */......？

所以我想我的问题是,为什么会发生这种情况以及实际发生了什么？我认为只要我可以分配所有相关的内存我应该没问题？对子程序的调用是否会复制参数,或只是指向它们的指针？如果答案是副本,那么我可以看到问题可能在哪里,如果是这样的话:关于如何绕过它的任何想法？

我试图解决的问题很大,但不以任何方式疯狂.标准的FE解算器可以处理比我现在更大的问题.我在Dell PowerEdge 1850上运行程序,操作系统是Microsoft Server 2008 R2 Enterprise.据systeminfo在cmd提示我的物理内存8GB和16GB几乎虚.据我所知,我的所有阵列和矩阵的总和不应超过100MB - 大约5.5M integer(4)和2.5M real(8)(根据我应该只有大约44MB,但让我们公平,并增加50MB的开销).

我使用与Microsoft Visual Studio 2008集成的英特尔Fortran编译器.

---

添加一些实际的源代码来澄清一点

! Update continuum state
CALL UpdateContinuumState(iTask,iArray,posc,dof,dof_k,nodedof,elm,&
                    bmtrx,detjac,w,mtrlprops,demtrx,dt,stress,strain,effstrain,&
                    effstress,aa,fi,errmsg)

是对例程的实际调用.大阵列是posc,bmtrx而且aa- 所有其他阵列至少要小一个数量级(如果不是更多).posc是INTEGER(4)和bmtrx而aa为REAL(8)

SUBROUTINE UpdateContinuumState(iTask,iArray,posc,dof,dof_k,nodedof,elm,bmtrx,&
                    detjac,w,mtrlprops,demtrx,dt,stress,strain,effstrain,&
                    effstress,aa,fi,errmsg)

    IMPLICIT NONE

    !I/O
    INTEGER(4) :: iTask, errmsg
    INTEGER(4) :: iArray(64)
    INTEGER(4),DIMENSION(iArray(15),iArray(15),iArray(5)) :: posc
    INTEGER(4),DIMENSION(iArray(22),iArray(21)+1) :: nodedof
    INTEGER(4),DIMENSION(iArray(29),iArray(3)+2) :: elm
    REAL(8),DIMENSION(iArray(14)) :: dof, dof_k
    REAL(8),DIMENSION(iArray(12)*iArray(17),iArray(15)*iArray(5)) :: bmtrx
    REAL(8),DIMENSION(iArray(5)*iArray(17)) :: detjac
    REAL(8),DIMENSION(iArray(17)) :: w
    REAL(8),DIMENSION(iArray(23),iArray(19)) :: mtrlprops
    REAL(8),DIMENSION(iArray(8),iArray(8),iArray(23)) :: demtrx
    REAL(8) :: dt
    REAL(8),DIMENSION(2,iArray(12)*iArray(17)*iArray(5)) :: stress
    REAL(8),DIMENSION(iArray(12)*iArray(17)*iArray(5)) :: strain
    REAL(8),DIMENSION(2,iArray(17)*iArray(5)) :: effstrain, effstress
    REAL(8),DIMENSION(iArray(25)) :: aa
    REAL(8),DIMENSION(iArray(14)) :: fi 

    !Locals
    INTEGER(4) :: i, e, mtrl, i1, i2, j1, j2, k1, k2, dim, planetype, elmnodes, &
        Nec, elmpnodes, Ndisp, Nstr, Ncomp, Ngpt, Ndofelm
    INTEGER(4),DIMENSION(iArray(15)) :: doflist
    REAL(8),DIMENSION(iArray(12)*iArray(17),iArray(15)) :: belm
    REAL(8),DIMENSION(iArray(17)) :: jelm
    REAL(8),DIMENSION(iArray(12)*iArray(17)*iArray(5)) :: dstrain
    REAL(8),DIMENSION(iArray(12)*iArray(17)) :: s
    REAL(8),DIMENSION(iArray(17)) :: ep, es, dep
    REAL(8),DIMENSION(iArray(15),iArray(15)) :: kelm
    REAL(8),DIMENSION(iArray(15)) :: felm

    dim       = iArray(1)
...

它在上一行之前失败了.

Answer 1

根据steabert的要求,我将在这里的评论中总结一下这里的对话更加明显,尽管MSB的答案已经正确解决了问题的核心问题.

在技​​术编程中,程序通常具有用于中间计算的大型本地数组,这发生了很多.局部变量通常存储在堆栈中,这通常(并且相当合理地)是整个系统内存的一小部分 - 通常为10MB左右.当局部变量大小超过堆栈大小时,您可以看到此处描述的症状 - 在调用相关子例程之后但在其第一个可执行语句之前发生堆栈溢出.

因此,当这个问题发生时,最好的办法是找到相关的大型局部变量,并决定做什么.在这种情况下,至少变量belm和dstrain变得非常大.

一旦找到变量,并且您已确认问题存在,那么有几个选项.正如MSB所指出的,如果你可以让你的阵列更小,那就是一个选择.或者,您可以使堆栈大小更大; 在linux下,完成了ulimit -s [newsize].但这确实只是推迟了问题,你必须在Windows机器上做一些不同的事情.

避免此问题的另一类方法不是将大数据放在堆栈上,而是放在内存的其余部分("堆")中.你可以通过赋予数组save属性(在C中static)来做到这一点; 这会将变量放在堆上,从而使值在调用之间保持不变.缺点是这可能会改变子程序的行为,并且意味着子程序不能递归使用,同样也是非线程安全的(如果你曾经处于多线程将同时进入例程的位置,它们每个人都会看到本地变量的相同副本,并可能会覆盖彼此的结果.好处是它简单易用 - 它应该可以在任何地方使用.但是,这只适用于固定大小的局部变量; 如果临时数组的大小取决于输入,则不能这样做(因为不再需要保存单个变量;每次调用过程时它的大小可能不同).

有特定于编译器的选项将所有数组(或所有大于某个给定大小的数组)放在堆而不是堆栈上; 我知道的每个Fortran编译器都有一个选项.对于在OPs帖子中使用的ifort,它-heap-arrays在linux或/heap-arrayswindows中.对于gfortran,这实际上可能是默认值.这有助于确保您知道发生了什么,但这意味着您必须为每个编译器提供不同的咒语以确保您的代码正常工作.

最后,您可以使有问题的数组可分配.分配的内存在堆上; 但是指向它们的变量在堆栈中,因此您可以获得两种方法的好处.此外,这是完全标准的fortran,因此非常便携.缺点是它需要更改代码.此外,分配过程可能需要花费大量时间; 因此,如果你打算召唤常规的数万次,你可能会注意到这会使事情略微减慢.(但是这种可能的性能回归很容易解决;如果你用相同大小的数组调用它数十亿次,你可以有一个可选的参数来传入预先分配的本地数组并使用它来代替,这样你只分配/解除分配一次).

每次分配/解除分配看起来像:

SUBROUTINE UpdateContinuumState(iTask,iArray,posc,dof,dof_k,nodedof,elm,bmtrx,&
                    detjac,w,mtrlprops,demtrx,dt,stress,strain,effstrain,&
                    effstress,aa,fi,errmsg)

    IMPLICIT NONE

    !...arguments.... 


    !Locals
    !...
    REAL(8),DIMENSION(:,:), allocatable :: belm
    REAL(8),DIMENSION(:), allocatable :: dstrain

    allocate(belm(iArray(12)*iArray(17),iArray(15))  
    allocate(dstrain(iArray(12)*iArray(17)*iArray(5))

    !... work

    deallocate(belm)
    deallocate(dstrain)

请注意,如果子例程执行了大量工作(例如,执行需要几秒钟),则一对分配/解除分配的开销应该是可以忽略的.如果没有,并且您希望避免开销,则使用preallocated worskpace的可选参数将类似于:

SUBROUTINE UpdateContinuumState(iTask,iArray,posc,dof,dof_k,nodedof,elm,bmtrx,&
                    detjac,w,mtrlprops,demtrx,dt,stress,strain,effstrain,&
                    effstress,aa,fi,errmsg,workbelm,workdstrain)

    IMPLICIT NONE

    !...arguments.... 
    real(8),dimension(:,:), optional, target :: workbelm
    real(8),dimension(:), optional, target :: workdstrain
    !Locals
    !...

    REAL(8),DIMENSION(:,:), pointer :: belm
    REAL(8),DIMENSION(:), pointer :: dstrain

    if (present(workbelm)) then
       belm => workbelm
    else
       allocate(belm(iArray(12)*iArray(17),iArray(15))
    endif
    if (present(workdstrain)) then
       dstrain => workdstrain
    else
       allocate(dstrain(iArray(12)*iArray(17)*iArray(5))
    endif

    !... work

    if (.not.(present(workbelm))) deallocate(belm)
    if (.not.(present(workdstrain))) deallocate(dstrain)