1

我想编写类似于 R 函数 seq() 的 Fortran 代码。例如：

x <- seq(0,1,0.1)

将给出向量

x <- c(0, 0.1, 0.2, ..., 1)

我将运行几个模拟，在这些模拟中序列的长度会改变，在 R 中这很容易做到，只需改变 seq() 中的第二个参数。我曾尝试在 Fortran 中使用动态数组和函数 ALLOCATE 来动态更改数组的大小。到目前为止，这还没有奏效并导致错误

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0  0x2B371ED7C7D7
#1  0x2B371ED7CDDE
#2  0x2B371F3B8FEF
#3  0x401BE9 in MAIN__ at test3D.f90:?
Segmentation fault (core dumped)

所以我想知道是否有一种简单的方法可以模仿 Fortran 中 R 函数 seq() 的行为。

如需进一步参考，请参阅下面的程序

program ffl
implicit none
integer, parameter           :: n = 2**12                  
integer                      :: m,j,l,o,num,r,posi         
real(kind=8), dimension(n)   :: results 
real(kind=8)                 :: dt,dk,dp, dtt, laenge, basal, periode,c      
real(kind=8), dimension(n,n) :: fitness, k_opt  
real(kind=8)                 :: t0,t1,t2,t3      
real(kind=8), dimension(:),allocatable    :: t   
real(kind=8), dimension(n)   :: k,p, tt1  
real(kind=8), dimension(6)   :: x_new, res, q0 
real(kind=8), dimension(6)   :: k1,k2,k3,k4    
real(kind=8)                 :: ts = 0.0    
real(kind=8)                 :: ks = 0.0, ke = 1.0  
real(kind=8)                 :: ps = 0.1, pe = 40.0  
real(kind=8)                 :: tts = 0.0, tte = 1.0  
real(kind=8), dimension(6)   :: u0,f1,f2,f3,u1    
external                     :: derivate 

! computing the vectors 
dk=(ke-ks)/real(n)    ! calculating resolution
dp=(pe-ps)/real(n)    ! calculating resolution
dtt=(tte-tts)/real(n) ! calculating resolution
k(1) = ks             ! first value for k = 0.0
p(1) = ps             ! first value for p = 0.001
tt1(1) = tts          ! first value for tts = 0.0

num = 10

do m = 1,n         
    k(m) = k(m-1)+dk ! setting the basal expression vector with resolution dt 
    tt1(m) = tt1(m-1)+dtt
end do

do m = 1,n
    p(m) = ps + 0.1
end do

do m = 1,n
    periode = p(m)

    do j = 1,n
    laenge = tt1(j)

        do l = 1,n
        basal = k(l)

            c = num * periode    ! calculating the length of the simulation
            dt=(c-ts)/real(n)    ! calculating time resolution
            r = 1
            t(1) = ts            ! setting first time value to t1 = 0

            allocate(t(1))       ! Initialize array dimension

            do while (ts + dt < c)
                t(r) = ts
                ts = ts + dt
                r = r + 1
                call resize_array
            end do

            ! initial conditions
            q0(1) = 0     ! x
            q0(2) = basal ! y
            q0(3) = 0     ! z
            q0(4) = 0     ! a
            q0(5) = 1     ! b
            q0(6) = 0     ! w 

            x_new = q0 ! set initial conditions
            ! Solving the model using a 4th order Runge-Kutta method
            do o = 1,n
                call derivate(basal,periode,laenge,t(l),x_new,k1)  

                t1 = t(o) + dt/2      
                f1 = x_new + (dt*k1)/2
                call derivate(basal,periode,laenge,t1,f1,k2)      

                t2 = t(o) + dt/2      
                f2 = x_new + (dt*k2)/2
                call derivate(basal,periode,laenge,t2,f2,k3)      

                t3 = t(o) + dt
                f3 = x_new + (dt*k3)/2
                call derivate(basal,periode,laenge,t3,f3,k4)      

                res = x_new + (dt*(k1+2*k2+2*k3+k4))/6
                if (res(2) < basal) then
                    res(2) = basal
                endif

                results(n) = res(6)

            end do
         fitness(j,l) = maxval(results)/c 
         end do
    write(*,*) fitness   
    !posi = maxloc(fitness(:,j)) 
    !k_opt(m,j) = k(posi)      ! inputting that value into the optimal k matrix
    end do
end do
!write(*,*) k_opt
!return k_opt

contains

! The subroutine increases the size of the array by 1
subroutine resize_array
real,dimension(:),allocatable :: tmp_arr
integer :: new

new = size(t) + 1

allocate(tmp_arr(new))
tmp_arr(1:new)=t
deallocate(t)

allocate(t(size(tmp_arr)))
t=tmp_arr

end subroutine resize_array   
end program ffl
4

2 回答 2

3

Fortran 2003 在分配可分配数组时具有（重新）分配，并且程序

program xgrid
implicit none
real, allocatable :: x(:)
integer           :: i,n
do n=5,10,5
   x = 0.1*[(i,i=0,n)]
   write (*,"('x =',100(1x,f0.1))") x
end do
end program xgrid

用 gfortran 4.8.0 编译，显示一个相当于 seq() 的 Fortran 单行，给出输出

x = .0 .1 .2 .3 .4 .5

x = .0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1.0

于 2014-04-23T13:21:37.103 回答
2

一个简单的实现，如果你真的想要这个函数，不想总是n手动计算。可能需要澄清一些上限。

print *,seq(1.,10.,0.1)

contains

function seq(from, to, step)
  real, allocatable :: seq(:)
  real, intent(in) :: from, to, step

  allocate(seq(0:int((to - from)/step)))
  do i = 0, size(seq)
    seq(i) = from + i * step
  end do
end function

end

关于您的程序，当您使用编译器具有的功能时，回溯会更有帮助。你的 resize_array 应该可能有tmp_arr(1:new-1)=t. move_alloc()子程序可以使它更短一点。

于 2014-04-23T13:39:12.143 回答