solve2.f90

运载火箭用拉瓦尔喷管计算（此程序参数保密

	SUBROUTINE DNETN(N,X,Y,EPS,FS,T,H,Aa,Bb,L,JS)
	DIMENSION X(N),Y(N),Aa(N,N),Bb(N),JS(N)
	DOUBLE PRECISION X,Y,Aa,Bb,T,H,AM,Z,BETA,D
	L=1000
10	CALL FS(X,Bb,N)
	AM=0.0
	DO 20 I=1,N
	  IF (ABS(Bb(I)).GT.AM) AM=ABS(Bb(I))
20	CONTINUE
	IF (AM.LT.EPS) RETURN
	L=L-1
	IF (L.EQ.0) THEN
	  WRITE(*,100)
	  RETURN
	END IF
100	FORMAT(1X,'  FAIL')
	DO 40 J=1,N
	  Z=X(J)
	  X(J)=X(J)+H
	  CALL FS(X,Y,N)
	  DO 30 I=1,N
30	  Aa(I,J)=Y(I)
	  X(J)=Z
40	CONTINUE
	CALL AGAUS(Aa,Bb,N,Y,K,JS)
	IF (K.EQ.0) THEN
	  L=-1
	  RETURN
	END IF
	BETA=1.0
	DO 50 I=1,N
50	BETA=BETA-Y(I)
	IF (ABS(BETA)+1.0.EQ.1.0) THEN
	  L=-2
	  WRITE(*,100)
	  RETURN
	END IF
	D=H/BETA
	DO 60 I=1,N
60	X(I)=X(I)-D*Y(I)
	H=T*H
	GOTO 10
	END


SUBROUTINE AGAUS(Aa,Bb,N,X,L,JS)
	DIMENSION Aa(N,N),X(N),Bb(N),JS(N)
	DOUBLE PRECISION Aa,Bb,X,T
	L=1
	DO 50 K=1,N-1
	  D=0.0
	  DO 210 I=K,N
	  DO 210 J=K,N
	    IF (ABS(Aa(I,J)).GT.D) THEN
	      D=ABS(Aa(I,J))
	      JS(K)=J
	      IS=I
	    END IF
210	  CONTINUE
	  IF (D+1.0.EQ.1.0) THEN
	    L=0
	  ELSE
	    IF (JS(K).NE.K) THEN
	      DO 220 I=1,N
	        T=Aa(I,K)
	        Aa(I,K)=Aa(I,JS(K))
	        Aa(I,JS(K))=T
220	      CONTINUE
	    END IF
	    IF (IS.NE.K) THEN
	      DO 230 J=K,N
	        T=Aa(K,J)
	        Aa(K,J)=Aa(IS,J)
	        Aa(IS,J)=T
230	      CONTINUE
	      T=Bb(K)
	      Bb(K)=Bb(IS)
	      Bb(IS)=T
	    END IF
	  END IF
	  IF (L.EQ.0) THEN
	    WRITE(*,100)
	    RETURN
	  END IF
	  DO 10 J=K+1,N
	    Aa(K,J)=Aa(K,J)/Aa(K,K)
10	  CONTINUE
	  Bb(K)=Bb(K)/Aa(K,K)
	  DO 30 I=K+1,N
	    DO 20 J=K+1,N
	      Aa(I,J)=Aa(I,J)-Aa(I,K)*Aa(K,J)
20	    CONTINUE
	    Bb(I)=Bb(I)-Aa(I,K)*Bb(K)
30	  CONTINUE
50	CONTINUE
	IF (ABS(Aa(N,N))+1.0.EQ.1.0) THEN
	  L=0
	  WRITE(*,100)
	  RETURN
	END IF
	X(N)=Bb(N)/Aa(N,N)
	DO 70 I=N-1,1,-1
	  T=0.0
	  DO 60 J=I+1,N
	    T=T+Aa(I,J)*X(J)
60	  CONTINUE
	  X(I)=Bb(I)-T
70	CONTINUE
100	FORMAT(1X,' FAIL ')
	JS(N)=N
	DO 150 K=N,1,-1
	  IF (JS(K).NE.K) THEN
	    T=X(K)
	    X(K)=X(JS(K))
	    X(JS(K))=T
	  END IF
150	CONTINUE
	RETURN
	END