2014-04-23 49 views
1

我想編寫一個像R函數seq()一樣工作的Fortran代碼。例如:動態大小的Fortran數組,與R函數一樣簡單seq()

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) 

所以我想知道是否有模仿R函數序列的Fortran語言行爲()的簡單方法。

爲了進一步參考見下文

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 

回答

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編譯,表演Fortran單行內容等同於seq(),給出輸出

x = .0 .1.2 .3 .4 .5

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

+0

謝謝,這正是我所需要的。 – StefanF

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 

關於你的程序,當你使用編譯器的fretures時,回溯會更有幫助。你的resize_array應該可能有tmp_arr(1:new-1)=t。子程序move_alloc()可能會縮短一點。