**************************************************************

* THIS EQUATION SOLVER IMPLEMENTS THE NEWTON METHOD

***************************************************************

LOGICAL IV[ALLOCATABLE](:)

INTEGER*2 INDX[ALLOCATABLE](:)

CHARACTER*3 SYMB[ALLOCATABLE](:),CH

REAL X[ALLOCATABLE](:),F[ALLOCATABLE](:),F1[ALLOCATABLE](:),

&D[ALLOCATABLE](:,:),VEC[ALLOCATABLE](:)

WRITE(*,*)' Give: (1) No. of Eqs, (2) No. of variables, (3) Input

&unit, (4) Output unit'

READ(*,*) N,NV,IN,IOUT

ALLOCATE (X(NV),F(N),F1(N),D(N,N),SYMB(NV),IV(NV),INDX(N),VEC(N))

IF(IN.EQ.0 .OR. IN.EQ.5) WRITE(IN,100) NV

100 FORMAT(' Give',I3,' Lines with:',/,3X,'(1) Symbol for var.(3 Ch)',

&/,3X,'(2) K or U for known, or unknown and',/,3X,'(3) Value')

J=0

DO 10 I=1,NV

READ(IN,110) SYMB(I),CH,X(I)

110 FORMAT(A3,1X,A1,1X,F10.0)

IF(CH.EQ.'U' .OR. CH.EQ.'u') THEN

IV(I)=.TRUE.

J=J+1

ELSE IF(CH.EQ.'K' .OR. CH.EQ.'k') THEN

IV(I)=.FALSE.

ELSE

WRITE(*,*)' Error in Input for variable',I

STOP

ENDIF

10 CONTINUE

IF(J.EQ.N) GO TO 12

WRITE(*,*)' You gave',N,' eqs., but',J,' unknowns'

STOP

12 NCT=0

15 CALL FUNCT(X,F)

J=0

DO 30 JJ=1,NV

IF(IV(JJ)) THEN

XX=X(JJ)

X(JJ)=1.005*X(JJ)

J=J+1

CALL FUNCT(X,F1)

DO 20 I=1,N

20 D(I,J)=(F1(I)-F(I))/(X(JJ)-XX)

X(JJ)=XX

ENDIF

30 CONTINUE

CALL SOLVEQ(N,N,D,F,INDX,VEC)

NCT=NCT+1

SUM=0.

J=0

DO 40 I=1,NV

IF(IV(I)) THEN

J=J+1

X(I)=X(I)-F(J)

SUM=SUM+ABS(F(J))

ENDIF

40 CONTINUE

WRITE(*,*)' NCT=',NCT,' SUM=',SUM

IF(NCT.LT.20 .AND. SUM.GT. .0001) GO TO 15

WRITE(IOUT,120)(I,SYMB(I),X(I),I=1,NV)

120 FORMAT(I5,1X,A3,' =',F10.3)

DEALLOCATE (X,F,F1,D,SYMB,IV,INDX,VEC)

END

***************************************************************

SUBROUTINE FUNCT(X,F)

REAL F(22),X(41)

DATA E,G2,P,AP/.005,64.4,8.6714174E-6,5.4541539E-3/

F(1)=X(41)-X(1)-X(3)-X(35) !Unknowns

F(2)=X(1)-X(2)-X(36) ! 1=Q2 15=H4 29=L1

F(3)=X(3)-X(4)-X(37) ! 2=Q3 16=H5 30=L2

F(4)=X(2)+X(4)-X(5)-X(38) ! 3=Q4 17=f1 31=L3

F(5)=X(12)-X(40)+X(6) ! 4=Q5 18=f2 32=L4

F(6)=X(13)-X(12)+X(7) ! 5=Q6 19=f3 33=L5

F(7)=X(14)-X(12)+X(9) ! 6=hf1 20=f4 34=L6

F(8)=X(15)-X(13)+X(8) ! 7=hf2 21=f5 35=QJ1

F(9)=X(16)-X(15)+X(11) ! 8=hf3 22=f6 36=QJ2

F(10)=X(7)+X(8)-X(10)-X(9) ! 9=hf4 Knowns 37=QJ3

RF1=1./SQRT(X(17)) ! 10=hf5 23=D1 38=QJ4

RF2=1./SQRT(X(18)) ! 11=hf6 24=D2 39=QJ5

RF3=1./SQRT(X(19)) ! 12=H1 25=D3 40=WS1

RF4=1./SQRT(X(20)) ! 13=H2 26=D4 41=Q1

RF5=1./SQRT(X(21)) ! 14=H3 27=D5

RF6=1./SQRT(X(22)) ! 28=D6

F(11)=X(6)-X(17)*X(29)*12./X(23)*(X(41)/(AP*X(23)**2))**2/G2

F(12)=X(7)-X(18)*X(30)*12./X(24)*(X(1)/(AP*X(24)**2))**2/G2

F(13)=X(8)-X(19)*X(31)*12./X(25)*(X(2)/(AP*X(25)**2))**2/G2

F(14)=X(9)-X(20)*X(32)*12./X(26)*(X(3)/(AP*X(26)**2))**2/G2

F(15)=X(10)-X(21)*X(33)*12./X(27)*(X(4)/(AP*X(27)**2))**2/G2

F(16)=X(11)-X(22)*X(34)*12./X(28)*(X(5)/(AP*X(28)**2))**2/G2

F(17)=RF1-1.14+2.*ALOG10(E/X(23)+P*X(23)*RF1/X(41))

F(18)=RF2-1.14+2.*ALOG10(E/X(24)+P*X(24)*RF2/X(1))

F(19)=RF3-1.14+2.*ALOG10(E/X(25)+P*X(25)*RF3/X(2))

F(20)=RF4-1.14+2.*ALOG10(E/X(26)+P*X(26)*RF4/X(3))

F(21)=RF5-1.14+2.*ALOG10(E/X(27)+P*X(27)*RF5/X(4))

F(22)=RF6-1.14+2.*ALOG10(E/X(28)+P*X(28)*RF6/X(5))

RETURN

END

**********************************************************

SUBROUTINE SOLVEQ(N,NP,A,B,INDX,VEC)

REAL A(NP,NP),B(NP),VEC(N)

INTEGER*2 INDX(NP),LL

DO 20 I=1,N

AMAX=0.

DO 10 J=1,N

IF(ABS(A(I,J)).GT.AMAX) AMAX=ABS(A(I,J))

10 CONTINUE

IF(AMAX.EQ.0.) STOP 'Singular matrix.'

20 VEC(I)=1./AMAX

DO 65 J=1,N

IF(J.EQ.1) GO TO 32

DO 30 I=1,J-1

TOL=A(I,J)

DO 25 K=1,I-1

25 TOL=TOL-A(I,K)*A(K,J)

A(I,J)=TOL

30 CONTINUE

32 AMAX=0.

DO 45 I=J,N

TOL=A(I,J)

IF(J.EQ.1) GO TO 40

DO 35 K=1,J-1

35 TOL=TOL-A(I,K)*A(K,J)

A(I,J)=TOL

40 TEMP=VEC(I)*ABS(TOL)

IF(TEMP.LT.AMAX) GO TO 45

MAX=I

AMAX=TEMP

45 CONTINUE

IF(J.EQ.MAX) GO TO 55

DO 50 K=1,N

TEMP=A(MAX,K)

A(MAX,K)=A(J,K)

50 A(J,K)=TEMP

VEC(MAX)=VEC(J)

55 INDX(J)=MAX

IF(J.EQ.N) GO TO 65

IF(A(J,J).EQ.0.) A(J,J)=1.E-10

TEMP=1./A(J,J)

DO 60 I=J+1,N

60 A(I,J)=A(I,J)*TEMP

65 CONTINUE

IF(A(N,N).EQ.0.) A(N,N)=1.E-10

II=0

DO 75 I=1,N

LL=INDX(I)

TOLL=B(LL)

B(LL)=B(I)

IF (II.NE.0)THEN

DO 70 J=II,I-1

70 TOLL=TOLL-A(I,J)*B(J)

ELSE IF(TOLL.NE.0.) THEN

II=I

ENDIF

75 B(I)=TOLL

DO 85 I=N,1,-1

TOLL=B(I)

IF(I.EQ.N) GO TO 85

DO 80 J=I+1,N

80 TOLL=TOLL-A(I,J)*B(J)

85 B(I)=TOLL/A(I,I)

RETURN

END