C***********************************************************************
C                                                                      *
C                       MULLER'S METHOD                                *
C                                                                      *
C***********************************************************************
C
C
C
C     TO FIND A SOLUTION TO F(X) = 0 GIVEN THREE APPROXIMATIONS X0, X1
C     AND X2:
C
C     INPUT X0,X1,X2; TOLERANCE TOL; MAXIMUM NUMBER OF ITERATIONS N0.
C
C     OUTPUT APPROXIMATE SOLUTION P OR MESSAGE OF FAILURE.
C
C     THIS IMPLEMENTATION ALLOWS FOR A SWITCH TO COMPLEX ARITHMETIC
C     THE COEFFICIENTS ARE STORED IN THE VECTOR A, SO THE DIMENSION
C     OF A MAY HAVE TO BE CHANGED
      COMPLEX Z(4),G(4),CH(3),CDEL1(2),CDEL,CB,CD,CE
      DIMENSION H(3),F(4),X(4),DEL1(2),A(10)
      LOGICAL OK
      INTEGER INP,OUP,FLAG
      CHARACTER NAME*14,NAME1*14,AA*1
      OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL')
      OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL')
      WRITE(6,*) 'This is the Muller Method.'
      OK = .FALSE.
10    IF ( .NOT. OK) THEN
         WRITE(6,*) 'Choice of input method: '
         WRITE(6,*) '1. Input entry by entry from keyboard '
         WRITE(6,*) '2. Input data from a text file '
         WRITE(6,*) 'Choose 1 or 2 please '
         WRITE(6,*) ' '
         READ(5,*)  FLAG
         IF( ( FLAG .GE. 1 ) .AND. ( FLAG .LE. 2 )) OK = .TRUE.
         GOTO 10
      ENDIF
      IF (FLAG .EQ. 1) THEN
         OK = .FALSE.
20       IF (.NOT. OK ) THEN
            WRITE(6,*) 'Input the degree n of the polynomial'
            WRITE(6,*) ' '
            READ(5,*) N
            IF (N .GT. 0 ) THEN
               OK = .TRUE.
               N=N+1
               WRITE(6,*) 'Input the coefficients of the'
               WRITE(6,*) 'polynomial in ascending order '
               WRITE(6,*) 'of exponent at the prompt '
               WRITE(6,*) ' '
               DO 30 I = 1, N
                  J=I-1
                  WRITE(6,*) 'Input A( ',J,' )'
                  WRITE(6,*) ' '
                  READ(5,*) A(I)
30             CONTINUE
            ELSE
               WRITE(6,*) 'Number must be a positive integer '
            ENDIF
            GOTO 20
         ENDIF
      ENDIF
      IF (FLAG .EQ. 2) THEN
         WRITE(6,*) 'Has a text file been created containing the '
         WRITE(6,*) 'coefficients in ascending order of exponent '
         WRITE(6,*) 'Enter Y or N - letter within quotes '
         WRITE(6,*) ' '
         READ(5,*)  AA
         IF (( AA .EQ. 'Y' ) .OR.( AA .EQ. 'y' )) THEN
            WRITE(6,*) 'Input the file name in the form - '
            WRITE(6,*) 'drive:name.ext  contained in quotes'
            WRITE(6,*) 'as example:   ''A:DATA.DTA'' '
            WRITE(6,*) ' '
            READ(5,*)  NAME
            INP = 4
            OPEN(UNIT=INP,FILE=NAME,ACCESS='SEQUENTIAL')
            OK = .FALSE.
40          IF (.NOT. OK) THEN
               WRITE(6,*) 'Input degree N '
               WRITE(6,*) ' '
               READ(5,*) N
               IF ( N .GT. 0) THEN
                  OK = .TRUE.
                  N=N+1
                  DO 50 I = 1, N
                     READ(4,*) A(I)
50                CONTINUE
               CLOSE(UNIT=4)
               ELSE
                  WRITE(6,*) 'Number must be a positive integer '
               ENDIF
               GOTO 40
            ENDIF
         ELSE
            WRITE(6,*) 'Please create the input file.'
            WRITE(6,*) 'The program will end so the input file can '
            WRITE(6,*) 'be created. '
            OK = .FALSE.
         ENDIF
      ENDIF
      IF(.NOT.OK) GOTO 400
11    OK = .FALSE.
12    IF (OK) GOTO 13
         WRITE(6,*) 'Input tolerance '
         WRITE(6,*) ' '
         READ(5,*) TOL
         IF (TOL.LE.0.0) THEN
            WRITE(6,*) 'Tolerance must be positive '
            WRITE(6,*) ' '
         ELSE
            OK = .TRUE.
         ENDIF
      GOTO 12
13    OK = .FALSE.
14    IF (OK) GOTO 15
         WRITE(6,*) 'Input maximum number of iterations '
         WRITE(6,*) '- no decimal point '
         WRITE(6,*) ' '
         READ(5,*) M
         IF ( M .LE. 0 ) THEN
           WRITE(6,*) 'Must be positive integer '
           WRITE(6,*) ' '
         ELSE
           OK = .TRUE.
         ENDIF
      GOTO 14
15    CONTINUE
      WRITE(6,*) 'Input the first of three starting values.'
      WRITE(6,*) ' '
      READ(5,*) X(1)
      WRITE(6,*) 'Input the second of three starting values.'
      WRITE(6,*) ' '
      READ(5,*) X(2)
      WRITE(6,*) 'Input the third of three starting values.'
      WRITE(6,*) ' '
      READ(5,*) X(3)
      WRITE(6,*) 'Select output destinations: '
      WRITE(6,*) '1. Screen '
      WRITE(6,*) '2. Text file '
      WRITE(6,*) 'Enter 1 or 2 '
      WRITE(6,*) ' '
      READ(5,*) FLAG
      IF ( FLAG .EQ. 2 ) THEN
         WRITE(6,*) 'Input the file name in the form - '
         WRITE(6,*) 'drive:name.ext'
         WRITE(6,*) 'with the name contained within quotes'
         WRITE(6,*) 'as example:   ''A:OUTPUT.DTA'' '
         WRITE(6,*) ' '
         READ(5,*) NAME1
         OUP = 3
         OPEN(UNIT=OUP,FILE=NAME1,STATUS='NEW')
      ELSE
         OUP = 6
      ENDIF
      WRITE(OUP,*) 'MULLER METHOD'
C     EVALUATE F USING HORNER'S METHOD AND STORE IN VECTOR F
      DO 60 I=1,3
         F(I)=A(N)
         DO 60 J=2,N
            K=N-J+1
60          F(I)=F(I)*X(I)+A(K)
      WRITE(OUP,5) (X(I),F(I),I=1,3)
C     VARIABLE ISW USED TO NOTE A SWITCH TO COMPLEX ARITHMETIC
C     ISW=0 MEANS REAL ARITHMETIC, AND ISW=1 MEANS COMPLEX ARITHMETIC
      ISW=0
C     STEP 1
      H(1)=X(2)-X(1)
      H(2)=X(3)-X(2)
      DEL1(1)=(F(2)-F(1))/H(1)
      DEL1(2)=(F(3)-F(2))/H(2)
      DEL=(DEL1(2)-DEL1(1))/(H(2)+H(1))
      I=3
C     STEP 2
100   IF (I.GT.M) GOTO 200
C        STEPS 3-7 FOR REAL ARITHMETIC
         IF(ISW.EQ.0) THEN
C           STEP 3
            B=DEL1(2)+H(2)*DEL
            D=B*B-4*F(3)*DEL
C           TEST TO SEE IF NEED COMPLEX ARITHMETIC
            IF(D.GE.0.0) THEN
C             REAL ARITHMETIC/ TEST TO SEE IF STRAIGHT LINE
               IF(ABS(DEL).LE.1.0E-20) THEN
C                 STRAIGHT LINE/ TEST TO SEE IF HORIZONTAL
                  IF(ABS(DEL1(2)).LE.1.0E-20) THEN
C                    HORIZONTAL LINE
                     WRITE(OUP,6)
                     GOTO 400
                  ENDIF
C                 STRAIGHT LINE BUT NOT HORIZONTAL
                  X(4)=(F(3)-DEL1(2)*X(3))/DEL1(2)
                  H(3)=X(4)-X(3)
               ELSE
C                 NOT STRAIGHT LINE
                  D=SQRT(D)
C                 STEP 4
                  E=B+D
                  IF(ABS(B-D).GT.ABS(E)) E=B-D
C                 STEP 5
                  H(3)=-2*F(3)/E
                  X(4)=X(3)+H(3)
               ENDIF
C              EVALUATE F(X(2))=F(4) BY HORNER'S METHOD
               F(4)=A(N)
               DO 70 J=2,N
                  K=N-J+1
70             F(4)=F(4)*X(4)+A(K)
               WRITE(OUP,8) I,X(4),F(4)
C              STEP 6
               IF(ABS(H(3)).LE.TOL) THEN
                  WRITE(OUP,9)
                  GOTO 400
               ENDIF
C              STEP 7
               DO 80 J=1,2
                  H(J)=H(J+1)
                  X(J)=X(J+1)
80                F(J)=F(J+1)
               X(3)=X(4)
               F(3)=F(4)
               DEL1(1)=DEL1(2)
               DEL1(2)=(F(3)-F(2))/H(2)
               DEL=(DEL1(2)-DEL1(1))/(H(2)+H(1))
            ELSE
C           SWITCH TO COMPLEX ARITHMETIC
               ISW=1
               DO 160 J=1,3
                  Z(J)=X(J)
160               G(J)=F(J)
               DO 170 J=1,2
                  CDEL1(J)=DEL1(J)
170               CH(J)=H(J)
               CDEL=DEL
            ENDIF
         ENDIF
C        STEPS 3-7 COMPLEX ARITHMETIC
         IF(ISW.EQ.1) THEN
C          TEST IF STRAIGHT LINE
            IF(CABS(CDEL).LE.1.0E-20) THEN
C             STRAIGHT LINE/ TEST IF HORIZONTAL
               IF(CABS(CDEL1(1)).LE.1.0E-20) THEN
C                 HORIZONTAL LINE
                  WRITE(OUP,6)
                  GOTO 400
               ENDIF
C              STRAIGHT LINE BUT NOT HORIZONTAL
               Z(4)=(G(4)-CDEL1(2)*Z(3))/CDEL1(2)
               CH(3)=Z(4)-Z(3)
            ELSE
C           NOT STRAIGHT LINE
C              STEP 3
               CB=CDEL1(2)+CH(2)*CDEL
               CD=CB*CB-4*G(3)*CDEL
               CD=CSQRT(CD)
C              STEP 4
               CE=CB+CD
               IF(CABS(CB-CD).GT.CABS(CE)) CE=CB-CD
C              STEP 5
               CH(3)=-2*G(3)/CE
               Z(4)=Z(3)+CH(3)
            ENDIF
C           EVALUATE G(X(I))=G(4) BY HORNER'S METHOD
            G(4)=A(N)
            DO 180 J=2,N
               K=N-J+1
180         G(4)=G(4)*Z(4)+A(K)
            WRITE(OUP,111) I,Z(4),G(4)
C           STEP 6
            IF(CABS(CH(3)).LE.TOL) THEN
               WRITE(OUP,9)
               GOTO 400
            ENDIF
C           STEP 7  CONTINUED
            DO 190 J=1,2
               CH(J)=CH(J+1)
               Z(J)=Z(J+1)
190            G(J)=G(J+1)
            Z(3)=Z(4)
            G(3)=G(4)
            CDEL1(1)=CDEL1(2)
            CDEL1(2)=(G(3)-G(2))/CH(2)
            CDEL=(CDEL1(2)-CDEL1(1))/(CH(2)+CH(1))
         ENDIF
C        STEP 7
         I=I+1
      GOTO 100
C     STEP 8
200   CONTINUE
      WRITE(OUP,121)
400   CLOSE(UNIT=5)
      CLOSE(UNIT=OUP)
      IF(OUP.NE.6) CLOSE(UNIT=6)
      STOP
3     FORMAT(8(1X,E15.8))
5     FORMAT(1X,E15.8,E15.8)
6     FORMAT(1X,'HORIZONTAL LINE')
8     FORMAT(1X,1X,I3,1X,E15.8,1X,E15.8)
9     FORMAT(1X,'PROCEDURE COMPLETED SUCCESSFULLY')
111   FORMAT(1X,I3,4(1X,E15.8))
121   FORMAT(1X,'FAILURE')
      END
