| 查看: 1329 | 回复: 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 |
» 猜你喜欢
时间戳又变了
已经有16人回复
科研孤儿太难了
已经有9人回复
时间戳变了,能看出什么问题?
已经有23人回复
只有每年这种时候来逛逛小木虫
已经有22人回复
今天系统多次维护,明天很可能放榜!
已经有10人回复
放榜前的不淡定
已经有12人回复
言之凿凿,明天官网会通知放榜?
已经有16人回复
重要消息,中午系统在维护
已经有12人回复
哪位老哥知道今年的国自然具体哪一天放榜?
已经有20人回复
今天基金会出结果吗?20260819
已经有15人回复
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
涂代洋
铁虫 (初入文坛)
- 应助: 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 :电子结构

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









回复此楼
40