[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode
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