| 查看: 1258 | 回复: 6 | ||
涂代洋铁虫 (初入文坛)
|
[求助]
fortran编程超出边界问题
|
|
|
program MAIN real*8 Hso,B,A,D,S,s6j,I,J,K,Y,X,P,N parameter(N=66) dimension B(N),A(N,N),D(N),S(N,N),Hso(N,N) CHARACTER*16 C(N) do i=1,N do j=1,N Hso(i,j)=0.d0 end do end do write(*,*)"Hso" do i=1,6 Hso(i,i)=((3.d0*4.d0*7.d0)**(1.0d0/2.0d0)) * *s6j(5.0,2.5,2.5,2.5,5.0,1.0) * *(39.d0**(1.0d0/2.0d0))*0.5d0 end do do i=7,14 Hso(i,i)=-((3.d0*4.d0*7.d0)**(1.0d0/2.0d0)) * *s6j(5.0,2.5,3.5,2.5,5.0,1.0) * *(39.d0**(1.0d0/2.0d0))*0.5d0 end do do i=15,24 Hso(i,i)=((3.d0*4.d0*7.d0)**(1.0d0/2.0d0)) * *s6j(5.0,2.5,4.5,2.5,5.0,1.0) * *(39.d0**(1.0d0/2.0d0))*0.5d0 end do do i=25,36 Hso(i,i)=-((3.d0*4.d0*7.d0)**(1.0d0/2.0d0)) * *s6j(5.0,2.5,5.5,2.5,5.0,1.0) * *(39.d0**(1.0d0/2.0d0))*0.5d0 end do do i=37,50 Hso(i,i)=((3.d0*4.d0*7.d0)**(1.0d0/2.0d0)) * *s6j(5.0,2.5,6.5,2.5,5.0,1.0) * *(39.d0**(1.0d0/2.0d0))*0.5d0 end do do i=51,66 Hso(i,i)=-((3.d0*4.d0*7.d0)**(1.0d0/2.0d0)) * *s6j(5.0,2.5,7.5,2.5,5.0,1.0) * *(39.d0**(1.0d0/2.0d0))*0.5d0 end dO do i=1, 66 do j=1, 66 A(i,j) = Hso(i,j) enddo enddo CALL JCB(N,A,S,1E-16) DO 30 I=1,N B(I)=A(I,I) 30 CONTINUE DO 40 J=1,N-1 P=J DO 50 I=J+1,N IF (B(I)<B(P))THEN P=I ENDIF 50 CONTINUE K=B(J) B(J)=B(P) B(P)=K 40 CONTINUE OPEN(2,FILE='E.DAT',STATUS='NEW') DO 60 I=1,N WRITE(2,200)I,B(I)-A(1,1) 60 CONTINUE OPEN(4,FILE='FS.DAT',STATUS='OLD') DO 110 I=1,N READ(4,300)C(I) 110 CONTINUE OPEN(7,FILE='VE.DAT',STATUS='NEW') DO 70 I=1,N DO 80 J=1,N IF(A(J,J).EQ.B(I))THEN WRITE(7,200)I,B(I)-A(1,1) DO 130 X=1,N D(X)=S(X,J) 130 CONTINUE DO 140 Y=1,N-1 P=Y DO 150 Z=Y+1,N IF (ABS(D(Z)).GT.ABS(D(P)))THEN P=Z ENDIF 150 CONTINUE K=D(Y) D(Y)=D(P) D(P)=K 140 CONTINUE P=1 DO 160 M=1,N DO 170 F=1,N IF(D(M).EQ.S(F,J).AND.D(M).NE.0)THEN WRITE(7,400)S(F,J),C(F) ELSEIF(D(M).EQ.S(M,J).AND.D(M).EQ.0)THEN WRITE(7,400)S(M,J),C(M) EXIT ELSEIF(D(M).EQ.0.AND.S(M,J).NE.0)THEN WRITE(7,400)D(M),C(P) P=P+1 EXIT ENDIF 170 CONTINUE 160 CONTINUE ENDIF 80 CONTINUE 70 CONTINUE 100 FORMAT(F25.16) 200 FORMAT(I4," ",F25.16) 300 FORMAT(A16) 400 FORMAT(F25.16," ",A16) CLOSE(4) CLOSE(2) CLOSE(7) stop end Function s6j(j1, j2, j3, l1, l2, l3) Real :: j1, j2, j3, l1, l2, l3, i, k,n, k1, k2,plus_p1,plus_sum real*8 s6j,p1 k1 = max(j1+j2+j3, j1+l2+l3, l1+j2+l3,l1+l2+l3) k2 = min(j1+j2+l1+l2, j2+j3+l2+l3, j3+j1+l3+l1) plus_sum = 0 Do k = k1, k2, 1 If (k-j1-j2-j3<0.or.abs(k-j1-j2-j3-int(k-j1-j2-j3))>0.0001.or. *k-j1-l2-l3<0.or.abs(k-j1-l2-l3-int(k-j1-l2-l3))>0.0001.or. *k-l1-j2-l3<0.or.abs(k-l1-j2-l3-int(k-l1-j2-l3))>0.0001.or. *k-l1-l2-j3<0.or.abs(k-l1-l2-j3-int(k-l1-l2-j3))>0.0001.or. *j1+j2+l1+l2-k<0.or.abs(j1+j2+l1+l2-k-int(j1+j2+l1+l2-k))>0.001.or. *j2+j3+l2+l3-k<0.or.abs(j2+j3+l2+l3-k-int(j2+j3+l2+l3-k))>0.001.or. *j3+j1+l3+l1-k<0.or.abs(j3+j1+l3+l1-k-int(j3+j1+l3+l1-k))>0.01)Then plus_p1 =0 Else plus_p1 = (((-1)**(k))*p1(k+1))/(p1(k-j1-j2-j3)* *p1(k-j1-l2-l3)*p1(k-l1-j2-l3)*p1(k-l1-l2-j3)* *p1(j1+j2+l1+l2-k)*p1(j2+j3+l2+l3-k)*p1(j3+j1+l3+l1-k)) plus_sum = plus_sum + plus_p1 End If End Do If (j1+j2-j3<0.or.abs(j1+j2-j3-int(j1+j2-j3))>0.0001.or. *j1-j2+j3<0.or.abs(j1-j2+j3-int(j1-j2+j3))>0.0001.or.-j1+j2+j3<0 *.or.abs(-j1+j2+j3-int(-j1+j2+j3))>0.0001.or.j1+j2+j3+1<0 *.or.abs(j1+j2+j3+1-int(j1+j2+j3+1))>0.or.j1+l2-l3<0.or. *abs(j1+l2-l3-int(j1+l2-l3))>0.0001.or.j1-l2+l3<0.or. *abs(j1-l2+l3-int(j1-l2+l3))>0.0001.or.-j1+l2+l3<0 *.or.abs(-j1+l2+l3-int(-j1+l2+l3))>0.0001.or.j1+l2+l3+1<0.or. *abs(j1+l2+l3+1-int(j1+l2+l3+1))>0.0001.or.l1+j2-l3<0.or. *abs(l1+j2-l3-int(l1+j2-l3))>0.0001.or.l1-j2+l3<0.or. *abs(l1-j2+l3-int(l1-j2+l3))>0.0001.or.-l1+j2+l3<0.or. *abs(-l1+j2+l3-int(-l1+j2+l3))>0.0001.or.l1+j2+l3+1<0.or. *abs(l1+j2+l3+1-int(l1+j2+l3+1))>0.0001.or.l1+l2-j3<0.or. *abs(l1+l2-j3-int(l1+l2-j3))>0.0001.or.l1-l2+j3<0.or. *abs(l1-l2+j3-int(l1-l2+j3))>0.0001.or.-l1+l2+j3<0.or. *abs(-l1+l2+j3-int(-l1+l2+j3))>0.0001.or.l1+l2+j3+1<0.or. *abs(l1+l2+j3+1-int(l1+l2+j3+1))>0.0001) THEN S6J=0 Else s6j=(((p1(j1+j2-j3)*p1(j1-j2+j3)*p1(-j1+j2+j3))/p1(j1+j2+j3+1)) * **(1.0/2.0))*(((p1(j1+l2-l3)*p1(j1-l2+l3)*p1(-j1+l2+l3)) * /p1(j1+l2+l3+1))**(1.0/2.0))*(((p1(l1+j2-l3)*p1(l1-j2+l3) * *p1(-l1+j2+l3))/p1(l1+j2+l3+1))**(1.0/2.0))*(((p1(l1+l2-j3) * *p1(l1-l2+j3)*p1(-l1+l2+j3))/p1(l1+l2+j3+1))**(1.0/2.0)) * *plus_sum End if End Function s6j Function p1(n) Real :: j1, j2, j3, l1, l2, l3, i, k,n, k1, k2 Real*8 p1 If (n==0) Then p1 = 1 Else p1 = 1 Do i = 1, n, 1 p1 = p1*i End Do End If End Function p1 SUBROUTINE JCB(N,A,S,EPS) real*8 A(N,N),S(N,N) DO 30 I=1,N DO 30 J=1,I IF(I-J) 20,10,20 10 S(I,J)=1.d0 GOTO 30 20 S(I,J)=0.d0 S(J,I)=0.d0 30 CONTINUE G=0.d0 DO 40 I=2,N I1=I-1 DO 40 J=1,I1 40 G=G+2.d0*A(I,J)*A(I,J) S1=SQRT(G) S2=EPS/FLOAT(N)*S1 S3=S1 L=0 50 S3=S3/FLOAT(N) 60 DO 130 IQ=2,N IQ1=IQ-1 DO 130 IP=1,IQ1 IF(ABS(A(IP,IQ)).LT.S3) GOTO 130 L=1 V1=A(IP,IP) V2=A(IP,IQ) V3=A(IQ,IQ) U=.5*(V1-V3) IF(U.EQ.0.) G=1.d0 IF(ABS(U).GE.1d-10) G=-SIGN(1.d0,U)*V2/SQRT(V2*V2+U*U) ST=G/SQRT(2.d0*(1.d0+SQRT(1.d0-G*G))) CT=SQRT(1.d0-ST*ST) DO 110 I=1,N G=A(I,IP)*CT-A(I,IQ)*ST !PRINT*,'^^^^^^^^^^^^^^^^^^^^^^^^^^^^' A(I,IQ)=A(I,IP)*ST+A(I,IQ)*CT A(I,IP)=G G=S(I,IP)*CT-S(I,IQ)*ST S(I,IQ)=S(I,IP)*ST+S(I,IQ)*CT 110 S(I,IP)=G DO 120 I=1,N A(IP,I)=A(I,IP) 120 A(IQ,I)=A(I,IQ) G=2.d0*V2*ST*CT A(IP,IP)=V1*CT*CT+V3*ST*ST-G A(IQ,IQ)=V1*ST*ST+V3*CT*CT+G A(IP,IQ)=(V1-V3)*ST*CT+V2*(CT*CT-ST*ST) A(IQ,IP)=A(IP,IQ) 130 CONTINUE IF(L-1) 150,140,150 140 L=0 GOTO 60 150 IF(S3.GT.S2) GOTO 50 RETURN END |
» 猜你喜欢
窗边初夏的小雨
已经有4人回复
2026年申博-电池方向
已经有11人回复
26年申博自荐-计算机视觉
已经有4人回复
导师各种操作恶心咋办
已经有8人回复
2026博士申请求助
已经有5人回复
研究生做的很差,你们会让毕业吗?
已经有11人回复
求碳排放博导;方向是LCA、生命周期可持续发展以及碳排放
已经有7人回复
2026博士或科研助理转27年博士
已经有7人回复
国自科送审了吗
已经有11人回复
博士招生
已经有5人回复
涂代洋
铁虫 (初入文坛)
- 应助: 0 (幼儿园)
- 金币: 5.3
- 散金: 20
- 帖子: 24
- 在线: 3.9小时
- 虫号: 2977137
- 注册: 2014-02-19
- 专业: 凝聚态物性 II :电子结构
5楼2016-07-13 10:21:55
argo
铁杆木虫 (著名写手)
上善若水
- 应助: 4 (幼儿园)
- 金币: 9672.8
- 散金: 200
- 红花: 34
- 帖子: 1940
- 在线: 486.4小时
- 虫号: 507533
- 注册: 2008-02-19
- 性别: GG
- 专业: 凝聚态物性 II :电子结构

2楼2016-07-12 12:53:35
argo
铁杆木虫 (著名写手)
上善若水
- 应助: 4 (幼儿园)
- 金币: 9672.8
- 散金: 200
- 红花: 34
- 帖子: 1940
- 在线: 486.4小时
- 虫号: 507533
- 注册: 2008-02-19
- 性别: GG
- 专业: 凝聚态物性 II :电子结构

3楼2016-07-12 13:11:27
涂代洋
铁虫 (初入文坛)
- 应助: 0 (幼儿园)
- 金币: 5.3
- 散金: 20
- 帖子: 24
- 在线: 3.9小时
- 虫号: 2977137
- 注册: 2014-02-19
- 专业: 凝聚态物性 II :电子结构
4楼2016-07-13 10:19:42
argo
铁杆木虫 (著名写手)
上善若水
- 应助: 4 (幼儿园)
- 金币: 9672.8
- 散金: 200
- 红花: 34
- 帖子: 1940
- 在线: 486.4小时
- 虫号: 507533
- 注册: 2008-02-19
- 性别: GG
- 专业: 凝聚态物性 II :电子结构

6楼2016-07-15 18:50:42
涂代洋
铁虫 (初入文坛)
- 应助: 0 (幼儿园)
- 金币: 5.3
- 散金: 20
- 帖子: 24
- 在线: 3.9小时
- 虫号: 2977137
- 注册: 2014-02-19
- 专业: 凝聚态物性 II :电子结构
7楼2016-07-16 11:08:28












回复此楼