我在这个Fortran 90代码中面临一个难以调试的问题-保留或注解掉write(*,*)
语句会改变代码的输出!下面是详细信息和代码(不能最小化更多,因为这可能不会导致相同的错误):
操作系统:Ubuntu 16.04
gfortran版本
GNU Fortran(Ubuntu 5.4.0-6ubuntu1~16.04.12)5.4.0 20160609
版权所有(C)2015自由软件基金会
FORTRAN 90代码名称:gillsp_notworking.f90
编译:gfortran -ffree-line-length-none gillsp_notworking.f90
下面是两个有问题的行(标记为!*error write):
W1:在第421行写入语句
W2:在第422行写入语句
代码:
program code
implicit none
real*8, parameter :: time=3600.d0, t0=0.d0
real*8, parameter :: er=1.d-6
real*8, parameter :: vton=600.d0
real*8, parameter :: dl=1.d0/370.d0
real*8,parameter :: kp1be = 11.6d0/vton, km1be = 1.4d0
real*8,parameter :: kp2be = 3.4d0/vton, km2be = 0.2d0
real*8,parameter :: kp3be = 2.9d0/vton, km3be = 5.4d0
real*8,parameter :: kp1pe = 0.d0/vton, km1pe = 0.d0
real*8,parameter :: kp2pe = 0.d0/vton, km2pe = 0.d0
real*8,parameter :: kp3pe = 0.d0/vton, km3pe = 0.d0
real*8,parameter :: kpdom = 16.d-5/vton
real*8,parameter :: kpcf = 13.d0/vton, kmcf = 0.7d0
real*8,parameter :: kmbcf = 4.4, kmpcf = 0.d0
real*8,parameter :: ktd = 0.003d0
real*8,parameter :: ktdpi = 0.3d0
real*8,parameter :: ksever = 0.01d0
real*8,parameter :: kcp0=12.8d0/vton, kcm0=2.d-4
real*8,parameter :: kcp1=0.3d0/vton, kcm1=6.d-3
real*8,parameter :: alpha = 20.d0
real*8,parameter :: rho = 0.3d0
real*8,parameter :: rhocf = 0.8d0
real*8,parameter :: rhocp = 0.d0
real*8,parameter :: V = 1.d0
integer, parameter :: Nmin=100,Npos=nint(60.d0/dl)
integer, parameter :: N=nint(rho*vton*V)
integer, parameter :: Ncf=nint(rhocf*vton*V)
integer, parameter :: Ncp=nint(rhocp*vton*V/1000.d0)
integer, parameter :: nstep=300
integer, parameter :: ndmax=3000
integer, parameter :: minsevL=10
integer, parameter :: plotmovie=1
real*8, parameter :: delt=time/nstep
real*8 :: t,r1,r2,pl,ph
real*8 :: tau
real*8 :: beta1,beta2,beta3,beta4,beta5,beta6,beta7,beta8,beta_sev
real*8 :: beta_cp,beta_cm
real*8 :: beta_sum
real*8 :: alpha0
real*8 :: rsev1,rsev2
real*8 :: pifactr
integer :: m, m1(-Nmin:Npos), m10(-Nmin:Npos)
integer :: mcf, domL(1:ndmax),domR(1:ndmax), ndom,cofcount,domLen
integer :: cap
integer :: domL0(1:ndmax),domR0(1:ndmax),drcount,domrem(1:ndmax)
integer :: pesever(1:10),besever(1:10),besvcount,pesvcount, sevindex, dumind
integer :: i,j,k,nk
integer :: l1(-Nmin:Npos)
integer :: l1bind, l1pind, l1bind0, l1pind0
integer :: gatp1, gadppi1, gadp1, glcf1, gzero
integer :: l1sum
integer :: k1,k2, dgcount, dscount, nmerge, domflag
integer :: nfile
integer :: dumdl
integer :: Bstate
integer :: nPicof
integer :: seed(12)
character*80::fname, Rid
write(*,*) 'Array size=', -Nmin, Npos
!------------------------------------------------------
seed= 9878771
CALL RANDOM_SEED (PUT=seed)
!------------------------------------------------------
t=0
nk=0
nfile=0
!----------- initialisation -----------
l1=0
m1=0
l1bind=0
l1pind=-1
m=N
cap=Ncp
Bstate = 1
pifactr = 1.d0
ndom=0
cofcount=0
domL=0
domR=0
nPicof=0
mcf=Ncf-sum(domR-domL)-ndom
!********************************************************************
!********************************************************************
!------------------- time loop ------------------------
tloop: do while (t.le.time)
!----------------- counting G-atp, G-adp, G-cofilin in filament --------------
gzero=count(m1<er)
gatp1=count(m1<1+er) - gzero
gadppi1=count(m1<2+er) - (gatp1+gzero)
gadp1=count(m1<3+er) - (gatp1+gzero+gadppi1)
glcf1=count(m1>3+er)
!----------------- current length ---------------------
l1sum = sum(l1)
!------------------ propensity sum ---------------------
!---- store monomer states, domain info etc -----
m10 = m1
l1bind0=l1bind
l1pind0=l1pind
domL0=domL
domR0=domR
!------------------- capping reaction ------------------
if (m10(l1bind-1).eq.4) then
if(Bstate.eq.1)then
beta_cp = kcp1*cap/V
beta_cm = 0.d0
elseif (Bstate.eq.2) then
beta_cp = 0.d0
beta_cm = kcm1
endif
else
if(Bstate.eq.1)then
beta_cp = kcp0*cap/V
beta_cm = 0.d0
elseif (Bstate.eq.2) then
beta_cp = 0.d0
beta_cm = kcm0
endif
endif
!-------------------- barbed end -----------------------
if(Bstate.eq.1)then ! Barbed-end state if
if(m10(l1bind-1).eq.1)then
beta1 = km1be
beta2 = kp1be*(m/V)
elseif (m10(l1bind-1).eq.2) then
beta1 = km2be
beta2 = kp2be*(m/V)
elseif (m10(l1bind-1).eq.3) then
beta1 = km3be
beta2 = kp3be*(m/V)
elseif (m10(l1bind-1).eq.4) then
beta1 = kmbcf
beta2 = 0.d0
endif
elseif (Bstate.eq.2) then
beta1 = 0.d0
beta2 = 0.d0
endif
!----- the first monomer incorporation------!
if(t.lt.30)then
if(m10(l1bind-1).eq.0)then
beta2 = kp1be*(m/V)
endif
endif
!-------------------- pointed end -----------------------
if(m10(l1pind+1).eq.1) then
beta3 = km1pe
beta4=kp1pe*(m/V)
elseif (m10(l1pind+1).eq.2) then
beta3 = km2pe
beta4=kp2pe*(m/V)
elseif (m10(l1pind+1).eq.3) then
beta3 = km3pe
beta4 = kp3pe*(m/V)
elseif (m10(l1pind+1).eq.4) then
beta3 = kmpcf
beta4 = 0.d0
endif
!------------------- bulk reactions ---------------------
!------------------ cofilin reactions ------------------
if(ndom.gt.0)then
nmerge=0
do k1=1,ndom
do k2=1,ndom
if(domL(k1)-1.eq.domR(k2)) nmerge=nmerge+1
enddo
enddo
endif
!---------------- Pi-cofilin interface ------------------
nPicof=0
do k1=l1pind0,l1bind0-1
if(abs(m10(k1)-m10(k1+1)).eq.2.and.m10(k1)+m10(k1+1).eq.6)then
nPicof=nPicof+1
endif
enddo
!------------------- bulk reactions ---------------------
beta5 = ktdpi*gatp1 + ktd*(gadppi1-nPicof) + alpha*ktd*nPicof
beta6 = kpdom*(mcf/V)*max(0.d0,(gadp1-2.d0*ndom-2.d0*nmerge))
!----- count domain growth probability -----
if(ndom.gt.0)then
dgcount=0
dscount=0
do k=1,ndom
dumdl=domR(k)-domL(k)-1
if(m10(domL(k)).eq.3) dgcount=dgcount+1
if(m10(domR(k)).eq.3) dgcount=dgcount+1
if(m10(domL(k)).le.3.and.dumdl.gt.2) dscount=dscount+1
if(m10(domR(k)).le.3.and.dumdl.gt.2) dscount=dscount+1
enddo
else
dgcount=0
dscount=0
endif
beta7 = kpcf*(mcf/V)*dgcount
beta8 = kmcf*dscount
!----- severing probability -----
if(ndom.gt.0)then
besvcount=0
pesvcount=0
besever=-1
pesever=-1
do k=1,ndom
domLen = domR(k)-domL(k)-1
if (domLen.ge.minsevL) then
if(m10(domL(k)).lt.4.and.m10(domL(k)).gt.0) then
pesvcount=pesvcount+1
pesever(pesvcount)=domL(k)
endif
if(m10(domR(k)).lt.4.and.m10(domR(k)).gt.0) then
besvcount=besvcount+1
besever(besvcount)=domR(k)
endif
else !------ domain length condition
pesvcount=0
besvcount=0
endif !------ domain length condition
enddo
else
pesvcount=0
besvcount=0
endif
beta_sev = ksever*(besvcount+pesvcount)
beta_sum = beta1+beta2+beta3+beta4+beta5+beta6+beta7 &
+beta8+beta_sev+beta_cp+beta_cm
alpha0 = beta_sum
!----------------- next reaction time: tau --------------
CALL RANDOM_NUMBER(r1)
if(r1.lt.1.d-8)then
CALL RANDOM_NUMBER(r1)
endif
tau = (1.d0/alpha0)*log(1.d0/r1)
!----------------------- all reactions ------------------------
CALL RANDOM_NUMBER(r2)
!---------------- growth and decay ------------------
!---------------- filament 1 ----------------
!---------------- barbed end ----------------
pl=0.d0
ph=beta1/alpha0
if (r2.ge.pl.and.r2.lt.ph.and.l1sum.gt.er) then
if (m1(l1bind-1).eq.4) then
! mcf = mcf+1 Pool conc. remains same
do k=1,ndom
if(domR(k).eq.l1bind) domR(k)=domR(k)-1
enddo
endif
l1bind = l1bind - 1
l1(l1bind) = 0
m1(l1bind) = 0
! m = m+1 Pool conc. remains same
Rid="BE-"
endif
pl=beta1/alpha0
ph=(beta1+beta2)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then
m1(l1bind) = 1
l1(l1bind) = 1
l1bind = l1bind + 1
! m = m-1 Pool conc. remains same
Rid="BE+"
endif
!---------------- pointed end ---------------
pl=(beta1+beta2)/alpha0
ph=(beta1+beta2+beta3)/alpha0
if (r2.ge.pl.and.r2.lt.ph.and.l1sum.gt.er) then
if (m1(l1pind+1).eq.4) then
! mcf = mcf+1 Pool conc. remains same
do k=1,ndom
if(domL(k).eq.l1pind) domL(k)=domL(k)+1
enddo
endif
l1pind = l1pind + 1
l1(l1pind) = 0
m1(l1pind) = 0
! m = m+1 Pool conc. remains same
Rid="PE-"
endif
pl=(beta1+beta2+beta3)/alpha0
ph=(beta1+beta2+beta3+beta4)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then
m1(l1pind) = 1
l1(l1pind) = 1
l1pind = l1pind - 1
! m = m-1 Pool conc. remains same
Rid="PE+"
endif
!---------- severing of filament -------------
pl=(beta1+beta2+beta3+beta4)/alpha0
ph=(beta1+beta2+beta3+beta4+beta_sev)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then
CALL RANDOM_NUMBER(rsev1)
CALL RANDOM_NUMBER(rsev2)
if(rsev1.le.0.2d0.and.besvcount.gt.0)then
dumind = 1 + nint(rsev2*(besvcount-1))
sevindex = besever(dumind)
else
if (pesvcount.gt.0) then
dumind = 1 + nint(rsev2*(pesvcount-1))
sevindex = pesever(dumind)
else
dumind = 1 + nint(rsev2*(besvcount-1))
sevindex = besever(dumind)
endif
endif
l1bind=sevindex ! keep olp pointed end
Bstate=1 ! set new barbed end free
! m=m+(l1bind0-sevindex) Pool conc. remains same
m1(sevindex:l1bind0)=0
l1(sevindex:l1bind0)=0
drcount=0
domrem=0
do j=sevindex,l1bind0 ! j-loop
! if(m10(j).eq.4) then Pool conc. remains same
! mcf=mcf+1
! endif
do k=1,ndom
if(domL(k).eq.j)then
drcount=drcount+1
domrem(drcount)=j
domL(k)=0
domR(k)=0
endif
enddo
enddo ! j-loop
Rid="SEVER"
! write(*,*)'sever', m10(sevindex-1) ,sevindex, l1bind0, besvcount, pesvcount !***** error write ****
! write(*,*)'sever', m10(sevindex-1) !***** error write ****
endif ! filament severing if
!---------------- capping ------------------
pl=ph
ph=(beta1+beta2+beta3+beta4+beta_sev+beta_cp)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then
! cap=cap-1 Pool conc. remains same
Bstate = 2
Rid="CAPPING +"
endif
!---------------- de-capping ------------------
pl=ph
ph=(beta1+beta2+beta3+beta4+beta_sev+beta_cp+beta_cm)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then
! cap=cap+1 Pool conc. remains same
Bstate = 1
Rid="CAPPING -"
endif
!---------- hydrolysis in filament -----------
pl=(beta1+beta2+beta3+beta4+beta_sev+beta_cp+beta_cm)/alpha0
do j=l1pind0,l1bind0
if(m10(j).eq.1)then
ph=pl+(ktdpi/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(j) = 2
Rid="HYD atp"
endif
pl=ph
elseif(m10(j).eq.2)then
if(m10(j-1).eq.4.or.m10(j+1).eq.4)then
pifactr=alpha
else
pifactr=1.d0
endif
ph=pl+(pifactr*ktd/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(j) = 3
Rid="Pi Release"
endif
pl=ph
elseif(m10(j).eq.3)then
domflag=0
if(ndom.gt.0)then ! ndom-if
do k=1,ndom
if(j.eq.domL(k)) then
domflag=1
endif
if(j.eq.domR(k)) then
domflag=1
endif
enddo
endif ! ndom-if
if (domflag.lt.1) then
ph=pl+((kpdom*mcf/V)/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(j) = 4
! mcf=mcf-1 Pool conc. remains same
ndom=ndom+1
domL(ndom)=j-1
domR(ndom)=j+1
Rid="DOM"
endif
pl=ph
endif
endif ! hydrolysis if
enddo ! j monomer loop
if(ndom.gt.0)then ! ndom if
do k=1,ndom
dumdl=domR(k)-domL(k)-1
!--------------------- domain growth --------------------
if (m10(domL(k)).eq.3) then ! Left boundary ADP-if
ph=pl+((kpcf*mcf/V)/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(domL(k))=4
! mcf=mcf-1 Pool conc. remains same
domL(k)=domL(k)-1
Rid="DOML+"
endif
pl=ph
endif
if (m10(domR(k)).eq.3) then ! Right boundary ADP-if
ph=pl+((kpcf*mcf/V)/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(domR(k))=4
! mcf=mcf-1 Pool conc. remains same
domR(k)=domR(k)+1
Rid="DOMR+"
endif
pl=ph
endif
! --------------------- domain shrinkage -------------------
if (m10(domL(k)).le.3.and.dumdl.gt.2) then
ph=pl+(kmcf/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(domL(k)+1)=3
! mcf=mcf+1 Pool conc. remains same
domL(k)=domL(k)+1
Rid="DOML-"
endif
pl=ph
endif
if (m10(domR(k)).le.3.and.dumdl.gt.2) then
ph=pl+(kmcf/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(domR(k)-1)=3
! mcf=mcf+1 Pool conc. remains same
domR(k)=domR(k)-1
Rid="DOMR-"
endif
pl=ph
endif
enddo ! k ndom loop
endif ! ndom if
!----------------- count bound cofilin --------------------
cofcount=0
do j=l1pind,l1bind
if (m1(j).eq.4) then
cofcount=cofcount+1
endif
enddo
!------------------- ERROR CHECKS ----------------------
if (tau.lt.0.d0) exit tloop
if (m1(l1bind).ne.0.or.m1(l1pind).ne.0) exit tloop
if (m.ne.N) exit tloop
if (mcf.ne.Ncf) exit tloop
if (cap.ne.Ncp) exit tloop
t = t + tau
enddo tloop !***** time loop *****
!------------------ time loop ended -----------------------
write(*,*) 'domains created', ndom
do k=1,ndom
write(*,*) k, 'size', domR(k)-domL(k)-1, 'L-R', domL(k), domR(k)
enddo
do k=l1pind,l1bind
write(101,*) k, m1(k)
enddo
if(t.ge.time) write(*,*) 'CODE ENDED WITHOUT ERROR'
write(*,*) 'CODE STOPPED ','t=',t,'id=',Rid
if (tau.lt.0.d0) write(*,*) 'ERROR: tau negative', tau, alpha0, beta1+beta2+beta3+beta4+beta5, beta6, beta7
if (m1(l1bind).ne.0.or.m1(l1pind).ne.0) write(*,*) 'ERROR: l1bind/l1pind not zero', m1(l1bind), m1(l1pind)
if (m.ne.N) write(*,*) 'ERROR: constant conc. fails', m, N
if (mcf.ne.Ncf) write(*,*) 'ERROR: constant conc. fails', mcf, Ncf
if (cap.ne.Ncp) write(*,*) 'ERROR: constant conc. fails', cap, Ncp
write(*,*) 'pe',l1pind,'be',l1bind, 'domains', ndom
stop
end program code
当W1和W2被注解掉时的O/P:
Array size= -100 22200
domains created 30
1 size 40 L-R 0 41
2 size 0 L-R -1 0
3 size -1 L-R 0 0
4 size -1 L-R 0 0
5 size -1 L-R 0 0
6 size -1 L-R 0 0
7 size 113 L-R 40 154
8 size -1 L-R 0 0
9 size -1 L-R 0 0
10 size -1 L-R 0 0
11 size -1 L-R 0 0
12 size 57 L-R 153 211
13 size -40 L-R 250 211
14 size -1 L-R 0 0
15 size -1 L-R 0 0
16 size -11 L-R 378 368
17 size -1 L-R 0 0
18 size -1 L-R 0 0
19 size -1 L-R 0 0
20 size -1 L-R 0 0
21 size -1 L-R 0 0
22 size -16 L-R 499 484
23 size -1 L-R 0 0
24 size -8 L-R 545 538
25 size -5 L-R 581 577
26 size -11 L-R 516 506
27 size -41 L-R 617 577
28 size -1 L-R 0 0
29 size -1 L-R 0 0
30 size -1 L-R 703 703
CODE ENDED WITHOUT ERROR
CODE STOPPED t= 3600.9739030733977 id=BE+
pe -1 be 239 domains 30
Note: The following floating-point exceptions are signalling: IEEE_DENORMAL
未注解W1时的O/P:
Array size= -100 22200
sever 4 144 550 2 1
sever 2 108 239 2 1
sever 2 207 475 1 1
sever 2 256 407 2 2
sever 2 224 579 2 2
sever 4 226 297 1 1
sever 2 337 562 2 2
sever 3 324 476 1 1
sever 4 376 611 1 1
sever 4 369 436 1 0
sever 4 615 943 1 1
sever 3 582 776 1 2
sever 4 409 734 1 1
sever 2 364 403 1 1
sever 4 512 833 3 2
sever 3 463 568 2 2
sever 4 450 493 1 1
sever 2 430 446 0 1
sever 3 865 1299 4 4
sever 4 812 1041 3 3
sever 4 547 813 1 1
sever 2 545 851 2 1
sever 2 573 765 1 1
domains created 31
1 size 122 L-R 0 123
2 size -1 L-R 0 0
3 size 0 L-R -1 0
4 size -1 L-R 0 0
5 size 75 L-R 122 198
6 size -1 L-R 0 0
7 size -1 L-R 0 0
8 size 39 L-R 197 237
9 size 48 L-R 236 285
10 size -1 L-R 0 0
11 size 3 L-R 316 320
12 size -1 L-R 0 0
13 size 0 L-R 356 357
14 size -1 L-R 0 0
15 size 36 L-R 319 356
16 size -1 L-R 0 0
17 size -1 L-R 0 0
18 size 139 L-R 357 497
19 size -1 L-R 0 0
20 size -1 L-R 0 0
21 size -1 L-R 0 0
22 size -1 L-R 0 0
23 size -1 L-R 0 0
24 size -1 L-R 0 0
25 size -1 L-R 0 0
26 size -1 L-R 0 0
27 size -1 L-R 0 0
28 size 47 L-R 517 565
29 size 21 L-R 496 518
30 size -1 L-R 0 0
31 size -1 L-R 0 0
CODE ENDED WITHOUT ERROR
CODE STOPPED t= 3600.0320511744094 id=HYD atp
pe -1 be 688 domains 31
Note: The following floating-point exceptions are signalling: IEEE_DENORMAL
没有一个结果是非物理的,不同的结果可以在Mac上用不同的编译器重现。我也用gfortran -Wall -Wno-tabs -g -fcheck=all -fbacktrace gillsp_notworking.f90
编译过,但没有发现任何严重的问题。除了不同的结果,代码有时也会在注解掉两个write语句的组合中冻结,这似乎取决于我之前运行的组合(疯狂!)-如果这对诊断有帮助的话。欢迎任何帮助或提示。提前谢谢!
------编辑1 -------
使用优化标志进行编译会抑制不同的输出,但在某些情况下会冻结(例如,当W1、W2未被注解掉时):gfortran -O gillsp_notworking.f90
并没有解决这个问题。
------编辑2 --------
根据评论的建议,我现在已经编译了调试选项和绑定检查选项:gfortran -Wall -Wno-tabs -Wextra -fbounds-check gillsp_notworking.f90
但这只是给了我一些关于未使用参数的警告。
根据使用random_number()
的建议,我修改了当前版本的代码和后续的不同行为。如果我保留两个write语句,它就会冻结。如果我使用-Og
标志(即gfortran -Og gillsp_notworking_seed.f90
)编译,当W1,W2都被注解掉时,结果就会改变!
1条答案
按热度按时间toe950271#
使用GFortran 11.3(Ubuntu 22.04中包含的默认版本),使用
-Og -Wall
编译,然后通过在初始化块中设置来修复may be used uninitialized
警告不管这两个write语句是否被注解掉,也不管优化级别如何,我都得到相同的输出(尝试了
-O0
、-O2
、-Og
、-O3
和-O3 -ffast-math
)。-fcheck=all
和-fsanitize=undefined
以及-fsanitize=address
均未发现任何问题(当然,这并不能保证问题不存在,但至少给出了一点指示)。