C  PCPRB8.FOR  27 February 1991
C
      PROGRAM PCPRB8
C
C  The Freudenstein-Roth function.
C
C  Reference
C
C  F Freudenstein, B Roth,
C  Numerical Solutions of Nonlinear Equations,
C  Journal of the Association for Computing Machinery,
C  Volume 10, 1963, Pages 550-556.
C
C  The function F(X) is of the form
C
C  FX(1) = X1 - X2**3 + 5*X2**2 -  2*X2 - 13 + 34*(X3-1)
C  FX(2) = X1 + X2**3 +   X2**2 - 14*X2 - 29 + 10*(X3-1)
C
C  Starting from the point (15,-2,0), the program is required to produce
C  solution points along the curve until it reaches a solution point
C  (*,*,1).  It also may be requested to look for limit points in the
C  first or third components.
C
C  The correct value of the solution at X3=1 is (5,4,1).
C
C  Limit points in the first variable occur at:
C
C  (14.28309, -1.741377,  0.2585779)
C  (61.66936,  1.983801, -0.6638797)
C
C  Limit points for the third variable occur at:
C
C  (20.48586, -0.8968053, 0.5875873)
C  (61.02031,  2.230139, -0.6863528)
C
      INTEGER   LIW
      INTEGER   LRW
      INTEGER   NVAR
C
      PARAMETER (NVAR=3)
      PARAMETER (LIW=NVAR+29)
      PARAMETER (LRW=29+(6+NVAR)*NVAR)
C
      EXTERNAL  DENSLV
      EXTERNAL  DFROTH
      EXTERNAL  FXROTH
      EXTERNAL  PITCON
C
      REAL      FPAR(1)
      INTEGER   I
      INTEGER   IERROR
      INTEGER   IPAR(1)
      INTEGER   ITRY
      INTEGER   IWORK(LIW)
      INTEGER   J
      INTEGER   LOUNIT
      CHARACTER NAME*12
      REAL      RWORK(LRW)
      REAL      XR(NVAR)
C
      LOUNIT=6
C
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'PCPRB8:'
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'PITCON sample program.'
      WRITE(LOUNIT,*)'Freudenstein-Roth function'
      WRITE(LOUNIT,*)'2 equations, 3 variables.'
C
C  Carry out continuation for each method of approximating the jacobian:
C
      ITRY=0
99    CONTINUE
C
C  Set work arrays to zero:
C
      DO 10 I=1,LIW
        IWORK(I)=0
10      CONTINUE
      DO 20 I=1,LRW
        RWORK(I)=0.0
20      CONTINUE
C
C  Set some entries of work arrays.
C
C  IWORK(1)=0 ; This is a startup
C  IWORK(2)=2 ; Use X(2) for initial parameter
C  IWORK(3)=0 ; Program may choose parameter index 
C  IWORK(4)=0 ; Update jacobian every newton step
C  IWORK(5)=3 ; Seek target values for X(3)
C  IWORK(6)=1 ; Seek limit points in X(1)
C  IWORK(7)=1 ; Control amount of output.
C  IWORK(8)=6 ; Output unit
C  IWORK(9)=* ; Jacobian choice.  Will vary from 0, 1, 2.
C               0=Use user's jacobian routine
C               1=Forward difference approximation
C               2=Central difference approximation
C
      IWORK(1)=0
      IWORK(2)=2
      IWORK(3)=0
      IWORK(4)=0
      IWORK(5)=3
      IWORK(6)=1
      IWORK(7)=1
      IWORK(8)=LOUNIT
      IWORK(9)=ITRY
C
C  RWORK(1)=0.00001; Absolute error tolerance
C  RWORK(2)=0.00001; Relative error tolerance
C  RWORK(3)=0.01   ; Minimum stepsize
C  RWORK(4)=20.0   ; Maximum stepsize
C  RWORK(5)=0.3    ; Starting stepsize
C  RWORK(6)=1.0    ; Starting direction
C  RWORK(7)=1.0    ; Target value (Seek solution with X(3)=1)
C
      RWORK(1)=0.00001
      RWORK(2)=0.00001
      RWORK(3)=0.01
      RWORK(4)=20.0
      RWORK(5)=0.3
      RWORK(6)=1.0
      RWORK(7)=1.0
C
C  Set starting point.
C
      XR(1)=15.0
      XR(2)=-2.0
      XR(3)=0.0
C
      WRITE(LOUNIT,*)' '
      IF(ITRY.EQ.0)THEN
        WRITE(LOUNIT,*)'Use user jacobian routine.'
      ELSEIF(ITRY.EQ.1)THEN
        WRITE(LOUNIT,*)'Approximate jacobian with forward difference.'
      ELSEIF(ITRY.EQ.2)THEN
        WRITE(LOUNIT,*)'Approximate jacobian with central difference.'
        ENDIF
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'Step  Type of Point     '//
     *'X(1)         X(2)         X(3)'
      WRITE(LOUNIT,*)' '
      I=0
      NAME='Start point  '
      WRITE(LOUNIT,'(1X,I3,2X,A12,2X,3G14.6)')I,NAME,(XR(J),J=1,NVAR)
C
C  Take another step.
C
      DO 30 I=1,30
C
        CALL PITCON(DFROTH,FPAR,FXROTH,IERROR,IPAR,IWORK,LIW,
     *  NVAR,RWORK,LRW,XR,DENSLV)
C
        IF(IWORK(1).EQ.1)THEN
          NAME='Corrected    '
        ELSEIF(IWORK(1).EQ.2)THEN
          NAME='Continuation '
        ELSEIF(IWORK(1).EQ.3)THEN
          NAME='Target point '
        ELSEIF(IWORK(1).EQ.4)THEN
          NAME='Limit point  '
          ENDIF
C
        WRITE(LOUNIT,'(1X,I3,2X,A12,2X,3G14.6)')I,NAME,(XR(J),J=1,NVAR)
C
        IF(IWORK(1).EQ.3)THEN
          WRITE(LOUNIT,*)' '
          WRITE(LOUNIT,*)'We have reached the point we wanted.'
          GO TO 40
          ENDIF
C
        IF(IERROR.NE.0)THEN
          WRITE(LOUNIT,*)' '
          WRITE(LOUNIT,*)'An error occurred.'
          WRITE(LOUNIT,*)'We will terminate the computation now.'
          GO TO 40
          ENDIF
C
30      CONTINUE
C
40    CONTINUE
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'Jacobians      ',IWORK(19)
      WRITE(LOUNIT,*)'Factorizations ',IWORK(20)
      WRITE(LOUNIT,*)'Solves         ',IWORK(21)
      WRITE(LOUNIT,*)'Functions      ',IWORK(22)
      ITRY=ITRY+1
      IF(ITRY.LE.2)GO TO 99
      STOP
      END
      SUBROUTINE FXROTH(NVAR,FPAR,IPAR,X,F,IERROR)
C
C  Evaluate the function at X.
C
C  ( X1 - ((X2-5.0)*X2+2.0)*X2 - 13.0 + 34.0*(X3-1.0)  )
C  ( X1 + ((X2+1.0)*X2-14.0)*X2 - 29.0 + 10.0*(X3-1.0) )
C
      INTEGER   NVAR
C
      REAL      F(*)
      REAL      FPAR(*)
      INTEGER   IERROR
      INTEGER   IPAR(*)
      REAL      X(NVAR)
C
      F(1)=X(1)
     *     -((X(2)-5.0)*X(2)+2.0)*X(2)
     *     -13.0
     *     +34.0*(X(3)-1.0)
C
      F(2)=X(1)
     *    +((X(2)+1.0)*X(2)-14.0)*X(2)
     *    -29.0
     *    +10.0*(X(3)-1.0)
C
      RETURN
      END
      SUBROUTINE DFROTH(NVAR,FPAR,IPAR,X,FJAC,IERROR)
C
C  Evaluate the Jacobian:
C
C  ( 1.0   (-3.0*X(2)+10.0)*X(2)-2.0      34.0   )
C  ( 1.0   (3.0*X(2)+2.0)*X(2)-14.0       10.0   )
C
      INTEGER   NVAR
C
      REAL      FJAC(NVAR,NVAR)
      REAL      FPAR(*)
      INTEGER   IERROR
      INTEGER   IPAR(*)
      REAL      X(NVAR)
C
      FJAC(1,1)=1.0
      FJAC(1,2)=(-3.0*X(2)+10.0)*X(2)-2.0
      FJAC(1,3)=34.0
C
      FJAC(2,1)=1.0
      FJAC(2,2)=(3.0*X(2)+2.0)*X(2)-14.0
      FJAC(2,3)=10.0
C
      RETURN
      END