| 查看: 1296 | 回复: 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 |
» 猜你喜欢
希望面上有个好结果
已经有5人回复
sci论文二审求助
已经有4人回复
今年也是没消息就是没中么
已经有12人回复
咸菜
已经有6人回复
有谁可曾问过你过的还好吗?
已经有5人回复
函评
已经有7人回复
以CTAB等为模板水热合成的产物,如何除去模板?
已经有6人回复
申博
已经有3人回复
前几天时间戳更新了
已经有14人回复
买卖文章的刷屏了!
已经有3人回复
涂代洋
铁虫 (初入文坛)
- 应助: 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











回复此楼