[Fortran] 纯文本查看 复制代码
EXTERNAL FS,F
DIMENSION JS(3),X(3)
DOUBLE PRECISION F,S,X
DATA JS/4,4,4/
N=3
CALL FGAUS(N,JS,X,FS,F,S)
WRITE(*,10) S
10 FORMAT(1X,'S=',D13.6)
END
SUBROUTINE FS(J,N,X,DN,UP)
DIMENSION X(N)
DOUBLE PRECISION X,DN,UP,Q
IF (J.EQ.1) THEN
DN=0.0
UP=1.0
ELSE IF (J.EQ.2) THEN
DN=0.0
UP=SQRT(1.0-X(1)*X(1))
ELSE IF (J.EQ.3) THEN
Q=X(1)*X(1)+X(2)*X(2)
DN=SQRT(Q)
UP=SQRT(2.0-Q)
END IF
RETURN
END
FUNCTION F(N,X)
DIMENSION X(N)
DOUBLE PRECISION F,X
F=X(3)*X(3)
RETURN
END
SUBROUTINE FGAUS(N,JS,X,FS,F,S)
DIMENSION JS(N),X(N)
DIMENSION T(5),C(5),D(2,11),CC(11),IS(2,11)
DOUBLE PRECISION X,F,S,T,C,D,CC,DN,UP,P
DATA T/-0.9061798459,-0.5384693101,0.0,
* 0.5384693101,0.9061798459/
DATA C/0.2369268851,0.4786286705,0.5688888889,
* 0.4786286705,0.2369268851/
M=1
D(1,N+1)=1.0
D(2,N+1)=1.0
10 DO 20 J=M,N
CALL FS(J,N,X,DN,UP)
D(1,J)=0.5*(UP-DN)/JS(J)
CC(J)=D(1,J)+DN
X(J)=D(1,J)*T(1)+CC(J)
D(2,J)=0.0
IS(1,J)=1
IS(2,J)=1
20 CONTINUE
J=N
30 K=IS(1,J)
IF (J.EQ.N) THEN
P=F(N,X)
ELSE
P=1.0
END IF
D(2,J)=D(2,J+1)*D(1,J+1)*P*C(K)+D(2,J)
IS(1,J)=IS(1,J)+1
IF (IS(1,J).GT.5) THEN
IF (IS(2,J).GE.JS(J)) THEN
J=J-1
IF (J.EQ.0) THEN
S=D(2,1)*D(1,1)
RETURN
END IF
GOTO 30
END IF
IS(2,J)=IS(2,J)+1
CC(J)=CC(J)+D(1,J)*2.0
IS(1,J)=1
END IF
K=IS(1,J)
X(J)=D(1,J)*T(K)+CC(J)
IF (J.EQ.N) GOTO 30
M=J+1
GOTO 10
END