C  PCPRB5.FOR  27 February 1991
C
      PROGRAM PCPRB5
C
C  Discretized two point boundary value problem with parameter LAMBDA:
C
C  Y''+LAMBDA*EXP(Y)=0
C  Y(0)=0 Y'(1)=0
C
C  Expecting limit point in LAMBDA at roughly LAMBDA=0.878.
C
C  This test illustrates the use of the full Newton, modified Newton,
C  and "cheap" Newton methods, which update the Jacobian on every iteration,
C  or only once per point, or only when convergence fails.
C
      INTEGER   LIW
      INTEGER   LRW
      INTEGER   NVAR
C
      PARAMETER (NVAR=22)
      PARAMETER (LRW=29+(NVAR+6)*NVAR)
      PARAMETER (LIW=NVAR+29)
C
      EXTERNAL  DENSLV
      EXTERNAL  FPFULL
      EXTERNAL  FXBEND
      EXTERNAL  PITCON
C
      REAL      FPAR(1)
      INTEGER   I
      INTEGER   IERROR
      INTEGER   IPAR(2)
      INTEGER   ITRY
      INTEGER   IWORK(LIW)
      INTEGER   J
      INTEGER   LOUNIT
      INTEGER   MODCON
      CHARACTER NAME*12
      INTEGER   NIT
      REAL      RWORK(LRW)
      REAL      XR(NVAR)
C
      LOUNIT=6
      ITRY=0
C
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'PCPRB5'
      WRITE(LOUNIT,*)'PITCON sample program.'
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'Two point BVP, limit point in Lambda'
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'This program will be run THREE times.'
C
19    CONTINUE
C
      ITRY=ITRY+1
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'This is run number ',ITRY
      NIT=0
C
      WRITE(LOUNIT,*)' '
      DO 4 I=1,LIW
        IWORK(I)=0
4       CONTINUE
      DO 5 I=1,LRW
        RWORK(I)=0.0
5       CONTINUE
      IERROR=0
C
C  Set input quantities
C
      IF(ITRY.EQ.1)MODCON=0
      IF(ITRY.EQ.2)MODCON=1
      IF(ITRY.EQ.3)MODCON=2
C
C  IWORK(1)=0      ; This is a startup
C  IWORK(2)=NVAR   ; Use X(NVAR) for initial parameter
C  IWORK(3)=0      ; Program may choose parameter index
C  IWORK(4)=MODCON ; Control frequency of jacobian update in newton iteration
C  IWORK(5)=NVAR   ; Seek target values for X(NVAR)
C  IWORK(6)=NVAR   ; Seek limit points in X(NVAR)
C  IWORK(7)=1      ; Control amount of output.
C  IWORK(8)=6      ; Output unit
C  IWORK(9)=0      ; Jacobian choice.
C
      IWORK(1)=0
      IWORK(2)=NVAR
      IWORK(3)=0
      IWORK(4)=MODCON
      IWORK(5)=NVAR
      IWORK(6)=NVAR
      IWORK(7)=1
      IWORK(8)=LOUNIT
      IWORK(9)=0
C
C  Jacobian upper/lower bandwidths:
C
      IPAR(1)=1
      IPAR(2)=1
C
C  RWORK(1)=0.0001; Absolute error tolerance
C  RWORK(2)=0.0001; Relative error tolerance
C  RWORK(3)=0.01   ; Minimum stepsize
C  RWORK(4)=0.25   ; Maximum stepsize
C  RWORK(5)=0.05   ; Starting stepsize
C  RWORK(6)=1.0    ; Starting direction
C  RWORK(7)=0.80   ; Target value (Seek solution with X(NVAR)=0.80
C
      RWORK(1)=0.0001
      RWORK(2)=0.0001
      RWORK(3)=0.01
      RWORK(4)=0.25
      RWORK(5)=0.05
      RWORK(6)=1.0
      RWORK(7)=0.80
C
C  Set starting point
C
      DO 55 I=1,NVAR
        XR(I)=0.0
55      CONTINUE
C
      IF(ITRY.EQ.1)THEN
        WRITE(LOUNIT,*)'Using full Newton method.'
        WRITE(LOUNIT,*)'Jacobian updated on every iteration.'
      ELSEIF(ITRY.EQ.2)THEN
        WRITE(LOUNIT,*)'Using modified Newton method.'
        WRITE(LOUNIT,*)'Jacobian held fixed while correcting a point.'
      ELSE
        WRITE(LOUNIT,*)'Using "cheap" Newton method.'
        WRITE(LOUNIT,*)'Jacobian held fixed as long as possible,'
        WRITE(LOUNIT,*)'perhaps over multiple points, until'
        WRITE(LOUNIT,*)'convergence fails.'
        ENDIF
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'Step  Type of Point     Lambda'
      WRITE(LOUNIT,*)' '
      I=0
      NAME='Start point  '
      WRITE(LOUNIT,'(1X,I3,2X,A12,2X,G14.6)')I,NAME,XR(NVAR)
C
      DO 10 I=1,50
        CALL PITCON(FPFULL,FPAR,FXBEND,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(NVAR)
C
        IF(IWORK(1).EQ.3)THEN
          NIT=NIT+1
          WRITE(LOUNIT,*)' '
          WRITE(LOUNIT,*)'Complete value for target point ',NIT
          WRITE(LOUNIT,*)' '
          WRITE(LOUNIT,'(1X,5G14.6)')(XR(J),J=1,NVAR)
          IF(NIT.GE.2)GO TO 20
          ENDIF
C
        IF(IERROR.NE.0)THEN
          WRITE(LOUNIT,*)'Iteration halted, IERROR=',IERROR
          GO TO 20
          ENDIF
C
10      CONTINUE
C
20    CONTINUE
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'Jacobians      ',IWORK(19)
      WRITE(LOUNIT,*)'Factorizations ',IWORK(20)
      WRITE(LOUNIT,*)'Solves         ',IWORK(21)
      WRITE(LOUNIT,*)'Functions      ',IWORK(22)
      IF(ITRY.LT.3)GO TO 19
      STOP
      END
      SUBROUTINE FXBEND(NVAR,FPAR,IPAR,X,FX,IERROR)
C
C***********************************************************************
C
C  Function evaluator for two point boundary value problem.
C
      INTRINSIC EXP
C
      INTEGER   NVAR
C
      REAL      EXP
      REAL      FPAR(1)
      REAL      FX(NVAR)
      REAL      HSQ
      INTEGER   I
      INTEGER   IERROR
      INTEGER   IPAR(2)
      REAL      X(NVAR)
C
      HSQ=NVAR-2
      HSQ=1.0/HSQ
      HSQ=HSQ*HSQ
      FX(1)=X(1)
      DO 10 I=2,NVAR-2
        FX(I)=X(I-1)-2.0*X(I)+X(I+1)+X(NVAR)*EXP(X(I))*HSQ
10      CONTINUE
      FX(NVAR-1)=X(NVAR-1)-X(NVAR-2)
      RETURN
      END
      SUBROUTINE FPFULL(NVAR,FPAR,IPAR,X,FPRIME,IERROR)
C
C***********************************************************************
C
C  Full storage jacobian of function
C
      INTRINSIC EXP
C
      INTEGER   NVAR
C
      REAL      EXP
      REAL      FPAR(1)
      REAL      FPRIME(NVAR,NVAR)
      REAL      HSQ
      INTEGER   I
      INTEGER   IERROR
      INTEGER   IPAR(2)
      REAL      X(NVAR)
C
      HSQ=NVAR-2
      HSQ=1.0/HSQ
      HSQ=HSQ*HSQ
      FPRIME(1,1)=1.0
      DO 10 I=2,NVAR-2
        FPRIME(I,I-1)=1.0
        FPRIME(I,I)=-2.0+X(NVAR)*EXP(X(I))*HSQ
        FPRIME(I,I+1)=1.0
        FPRIME(I,NVAR)=EXP(X(I))*HSQ
10      CONTINUE
      FPRIME(NVAR-1,NVAR-2)=-1.0
      FPRIME(NVAR-1,NVAR-1)=1.0
      RETURN
      END