C  PCPRB6.FOR  27 February 1991
C
      PROGRAM PCPRB6
C
C  Demonstrate use of IWORK(1)=-1 to check jacobian.
C
C  The problem solved is the same as problem number 1, the
C  Freudenstein-Roth function.
C
C  On the second try of the problem, the Jacobian is "polluted".
C  The Jacobian checker points this out.  Note, however, that the
C  calculation converges, all the way to the target point, including
C  limit points!  However, the number of function and jacobian evaluations
C  shoots up.  This illustrates the fact that an error in the jacobian
C  may not show up because of nonconvergence, but only because of
C  "slow" convergence!
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   ISAVE
      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,*)'PCPRB6'
      WRITE(LOUNIT,*)'PITCON sample program.'
      WRITE(LOUNIT,*)'Freudenstein-Roth function.'
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'Demonstrate use of the Jacobian checker'
      WRITE(LOUNIT,*)'on the third step.'
C
      ITRY=0
99    CONTINUE
C
C  Use IPAR(1) as a signal to Jacobian routine.  If IPAR(1)=1,
C  then use a 'bad' jacobian.
C
      ITRY=ITRY+1
      IF(ITRY.EQ.1)THEN
        IPAR(1)=0
      ELSE
        IPAR(1)=1
        ENDIF
C
      DO 10 I=1,LIW
        IWORK(I)=0
10      CONTINUE
      DO 20 I=1,LRW
        RWORK(I)=0.0
20      CONTINUE
C
C  IWORK(1)=0 ; This is a startup
C  IWORK(2)=2 ; Use index 2 for first parameter
C  IWORK(3)=0 ; Program may choose index
C  IWORK(4)=0 ; Update jacobian every newton step
C  IWORK(5)=3 ; Seek target values for index 3
C  IWORK(6)=1 ; Seek limit points in index 1
C  IWORK(7)=1 ; small amount of output
C  IWORK(8)=6 ; Output unit to use
C  IWORK(9)=0 ; Use user's jacobian routine
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)=0
C
C  RWORK(1)=0.00001; Absolute error tolerance
C  RWORK(2)=0.0001 ; 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.0001
      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.1)THEN
        WRITE(LOUNIT,*)'Use correct jacobian'
      ELSE
        WRITE(LOUNIT,*)'Repeat the run, but this time, use'
        WRITE(LOUNIT,*)'incorrect jacobian!'
        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
C  Check jacobian on third step
C
        IF(I.EQ.3)THEN
          WRITE(LOUNIT,*)' '
          WRITE(LOUNIT,*)'Turning on Jacobian check option!'
          WRITE(LOUNIT,*)' '
          ISAVE=IWORK(1)
          IWORK(1)=-1
          ENDIF
C
        CALL PITCON(DFROTH,FPAR,FXROTH,IERROR,IPAR,IWORK,LIW,
     *  NVAR,RWORK,LRW,XR,DENSLV)
C
C  After Jacobian check, restore value of IWORK(1)
C
        IF(I.EQ.3)THEN
          IWORK(1)=ISAVE
          ENDIF
        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
        IF(IERROR.NE.0)THEN
          WRITE(LOUNIT,*)'Iteration halted, IERROR=',IERROR
          GO TO 40
          ENDIF
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)
C
      IF(ITRY.LT.2)GO TO 99
      STOP
      END
      SUBROUTINE FXROTH(NVAR,FPAR,IPAR,X,F,IERROR)
C
      INTEGER   NVAR
C
      REAL      F(NVAR)
      REAL      FPAR(1)
      INTEGER   IERROR
      INTEGER   IPAR(1)
      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)
      RETURN
      END
      SUBROUTINE DFROTH(NVAR,FPAR,IPAR,X,FJAC,IERROR)
C
      INTEGER   NVAR
C
      REAL      FJAC(NVAR,NVAR)
      REAL      FPAR(1)
      INTEGER   IERROR
      INTEGER   IPAR(1)
      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
C  If IPAR(1)=1, put incorrect values in jacobian
C
      IF(IPAR(1).EQ.1)THEN
        FJAC(1,1)=1.2
        FJAC(1,2)=1.1*FJAC(1,2)
        FJAC(2,2)=0.9*FJAC(2,2)
        ENDIF
      RETURN
      END