fortran - Fortran 90中的堆栈溢出

Question

我已经在Fortran 90中编写了一个相当大的程序。它已经运行了一段时间了，但是今天我尝试将其提高一个档次并增加问题的大小(这是研究非标准的有限元求解器，如果那样的话)。可以帮助任何人...)现在，我收到“堆栈溢出”错误消息，程序自然终止，而没有提供任何有用的帮助。

该程序首先设置所有相关的数组和矩阵，然后完成此操作，并将与此相关的几行统计信息打印到日志文件中。即使出现了新的更大的问题，它也可以正常工作(尽管有点慢)，但是随着“数字紧缩”的进行，它就失败了。

令我感到困惑的是，那时所有的东西都已经分配了(并且可以正常工作而没有错误)。我不能完全确定堆栈是什么(维基百科和此处的几个步骤没有做太多事情，因为我对计算机的“幕后工作”只有一个非常基本的了解)。

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

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 */ ...？

所以我想我的问题是，为什么会发生这种情况以及实际上发生了什么？我以为只要分配所有相关的内存就可以了吗？对子例程的调用会复制参数，还是只是指向它们的指针？如果答案是副本，那么我可以看到问题可能在哪里，如果可以，那么如何解决该问题的任何想法？

我试图解决的问题很大，但无论如何都不是疯子。标准的有限元求解器可以处理比我当前更大的问题。我在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集成的Intel 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的要求，我将在这里的评论中总结一下对话，尽管M.S.B.的答案已经很贴切地解决了这个问题，但在这里，它看起来更为明显。

在技​​术编程中，过程通常具有较大的局部数组用于中间计算，这种情况经常发生。局部变量通常存储在堆栈中，通常(相当合理地)占整个系统内存的一小部分-通常约为10MB。当局部变量大小超过堆栈大小时，您会确切地看到此处描述的症状-在调用相关子例程之后但在其第一个可执行语句之前，发生堆栈溢出。

因此，当发生此问题时，最好的办法是找到相关的大型局部变量，然后决定要做什么。在这种情况下，至少变量belm和dstrain变得相当大。

找到变量后，您已经确认是问题所在，这里有一些选择。正如MSB指出的那样，如果可以使阵列更小，那是一种选择。另外，您可以增加堆栈大小；在linux下，这是通过ulimit -s [newsize]完成的。但是，这实际上只是推迟了该问题，您必须在Windows机器上做一些不同的事情。

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

有特定于编译器的选项，这些选项将所有数组(或所有大于某个给定大小的数组)放到堆而不是堆栈上。我知道的每个Fortran编译器都有一个选项。对于OPort中使用的ifort，在Linux中是-heap-arrays，在Windows中是/heap-arrays。对于gfortran，这实际上可能是默认设置。这有助于确保您知道发生了什么，但这意味着每个编译器都必须具有不同的要求，以确保您的代码能够正常工作。

最后，您可以使有问题的数组可分配。分配的内存在堆上；但是指向它们的变量在堆栈中，因此您可以从这两种方法中受益。而且，这是完全标准的fortran，因此完全可移植。缺点是需要更改代码。同样，分配过程可能要花费很短的时间。因此，如果您要调用无数次例行 session ，您可能会注意到这会稍微减慢速度。 (但是，这种可能的性能回归很容易修复；如果要使用相同大小的数组将其称为数十亿次，则可以有一个可选参数来传递预分配的本地数组，并使用该参数，以便您只能分配/取消分配一次)。

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

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)

请注意，如果子例程进行了大量工作(例如，需要花费几秒钟来执行)，那么来自几个分配/取消分配的开销应该可以忽略不计。如果不是这样，并且您想避免开销，那么为预分配的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)