Fortran排序法数值配方

Question

我使用C来调用Fortran，我的fortran正在调用sort()方法

 *-----------------------------------------------------------------------
 * SUBROUTINE sort(A,n)
 * Subroutine de la librairie "Numerical Recipes"
 * (C) Copr. 1986-92 Numerical Recipes Software
 *-----------------------------------------------------------------------
      SUBROUTINE sort(arr,n)
      INTEGER n,M,NSTACK
      REAL arr(n)
      PARAMETER (M=7,NSTACK=50)
      INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
      REAL a,temp
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
        do 12 j=l+1,ir
          a=arr(j)
          do 11 i=j-1,1,-1
            if(arr(i).le.a)goto 2
            arr(i+1)=arr(i)
11        continue
          i=0
2         arr(i+1)=a
12      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        temp=arr(k)
        arr(k)=arr(l+1)
        arr(l+1)=temp
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          arr(l+1)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l).gt.arr(ir))then
          temp=arr(l)
          arr(l)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l+1).gt.arr(l))then
          temp=arr(l+1)
          arr(l+1)=arr(l)
          arr(l)=temp
        endif
        i=l+1
        j=ir
        a=arr(l)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        arr(i)=arr(j)
        arr(j)=temp
        goto 3
5       arr(l)=arr(j)
        arr(j)=a
        jstack=jstack+2
        if(jstack.gt.NSTACK)pause 'NSTACK too small in sort'
        if(ir-i+1.ge.j-l)then
         istack(jstack)=ir
         istack(jstack-1)=i
         ir=j-1
       else
         istack(jstack)=j-1
         istack(jstack-1)=l
         l=i
       endif
     endif
     goto 1
     END

如果我多次调用排序方法，那么这个方法就有一个分段错误：

这是遗留代码，但我信任它，因为它来自数字配方。但我对一些事情很怀疑，尤其是这句话：

if(jstack.gt.NSTACK)pause 'NSTACK too small in sort'

如果我在这种情况下，我的节目会暂停？排序方法怎么可能做到这一点？

如果这一行可疑，我怎么能相信整个代码呢？

有人知道这类子程序有什么问题吗？有人知道用fortran进行排序的另一种方法吗？因为我可以用另一种方法代替这种方法，但是我是fortran的新手，我不能再写另一种方法。

我补充说，如果我在单线程中运行这个方法没有问题，但是如果我在多线程环境中运行它，问题就在这里。很抱歉，当我写我的问题时没有提到，但写完之后我看到了这个。

调试信息

current thread: t@41
  [1] __lwp_kill(0x0, 0x6, 0x0, 0x6, 0xffbffeff, 0x0), at 0xff2caa58
  [2] raise(0x6, 0x0, 0xff342f18, 0xff2aa378, 0xffffffff, 0x6), at 0xff265a5c
  [3] abort(0x7400, 0x1, 0x0, 0xfcb78, 0xff3413d8, 0x0), at 0xff24194c
  [4] os::abort(0x1, 0x0, 0xff011084, 0xfefdc000, 0x7d94, 0x7c00), at 0xfee7d3cc
  [5] VMError::report_and_die(0x0, 0xff038640, 0xff031ff4, 0x1, 0xfee81b94, 0xff031ff4), at 0xfef0cd58
  [6] JVM_handle_solaris_signal(0xb, 0xacffefe0, 0xacffed28, 0x8000, 0xff030fa0, 0x2013d8), at 0xfea73d48
  [7] __sighndlr(0xb, 0xacffefe0, 0xacffed28, 0xfea7325c, 0x0, 0x1), at 0xff2c6e78
  ---- called from signal handler with signal 11 (SIGSEGV) ------
  [8] sort_(0xfe2b1350, 0xfe2b135c, 0xfe2b1000, 0x1c00, 0x443bfc7b, 0xfe292484), at 0xfe27e498
  [9] mediane_(0xa9c1624c, 0xacfff2ac, 0xa9c16060, 0xa9c05c34, 0x0, 0x19), at 0xfe27a38c


(dbx) frame 8
0xfe27e498: sort_+0x01d8:       ld       [%l4 + %l1], %f4
(dbx) disassemble
0xff2caa58: __lwp_kill+0x0008:  bcc,a,pt  %icc,__lwp_kill+0x18  ! 0xff2caa68
0xff2caa5c: __lwp_kill+0x000c:  clr      %o0
0xff2caa60: __lwp_kill+0x0010:  cmp      %o0, 91
0xff2caa64: __lwp_kill+0x0014:  move     %icc,0x4, %o0
0xff2caa68: __lwp_kill+0x0018:  retl
0xff2caa6c: __lwp_kill+0x001c:  nop
0xff2caa70: __lwp_self       :  mov      164, %g1
0xff2caa74: __lwp_self+0x0004:  ta       %icc,0x00000008
0xff2caa78: __lwp_self+0x0008:  retl
0xff2caa7c: __lwp_self+0x000c:  nop

在solaris中的m中有dbx=> gdb

我试着检查入口，但是我能输入什么来获得有趣的信息呢？

在将-g选项添加到f90编译器后，在dbx中，我可以看到值或var，并看到结果：

t@88 (l@88) terminated by signal ABRT (Abort)
0xff2caa58: __lwp_kill+0x0008:  bcc,a,pt  %icc,__lwp_kill+0x18  ! 0xff2caa68
Current function is sort
  578           temp=arr(k)
(dbx) print n
n = 19
(dbx) print arr
arr =
    (1)    725.0666
    (2)    741.5034
    (3)    730.8196
    (4)    754.3707
    (5)    741.718
    (6)    741.718
    (7)    741.8914
    (8)    745.9141
    (9)    744.6705
    (10)    741.718
    (11)    745.8358
    (12)    743.3788
    (13)    746.2706
    (14)    746.2706
    (15)    750.1498
    (16)    754.3707
    (17)    754.3707
    (18)    754.3707
    (19)    748.2084
(dbx) print istack
istack =
    (1)    7
    (2)    12
    (3)    17
    (4)    18
    (5)    8
    (6)    9
    (7)    1
    (8)    4
    (9)    0
    (10)    0
    (11)    0
    (12)    0
    (13)    0
    (14)    0
    (15)    0
    (16)    0
    (17)    0
    (18)    0
    (19)    0
    (20)    0
    (21)    0
    (22)    0
    (23)    0
    (24)    0
    (25)    0
    (26)    0
    (27)    0
    (28)    0
    (29)    0
    (30)    0
    (31)    0
    (32)    0
    (33)    0
    (34)    0
    (35)    0
    (36)    0
    (37)    0
    (38)    0
    (39)    0
    (40)    0
    (41)    0
    (42)    0
    (43)    0
    (44)    0
    (45)    0
    (46)    0
    (47)    0
    (48)    0
    (49)    0
    (50)    0
(dbx) print jstack
jstack = -31648
(dbx)

它可能有一个-31648 val！istack只有50个元素，而istack(j堆栈)返回给我一个abd值！怎么可能呢？:)预先感谢

Answer 1

太大以至于不能发表评论：

不幸的是，堆栈跟踪没有显示出确切的失败线。错误发生在您所指示的线上的证据是什么？我可以想象在调用运行时系统时可能会出现错误，所以您应该尝试将PAUSE语句更改为write、read *,或类似的东西。

我在gfortran中用你的子程序做了一些测试：

  parameter (n = 100000)
  dimension b(n)

  do i=1,1000
    call random_number(b)
    call sort(b,n)
  end do
end

具有不同的数组大小和循环边界。这使用许多不同的输入调用sort。我启用了所有的检查和消毒，没有遇到任何一个问题。

编辑：

它也适用于OpenMP：

  parameter (n = 100000)
  real,allocatable :: b(:)

  !$omp parallel private(b)
  allocate(b(n))
  !$omp do
  do i=1,1000
    call random_number(b)
    call sort(b,n)
  end do
  !$omp end do
  !$omp end parallel
end