C  PCPRB3.FOR  27 February 1991
C
      PROGRAM PCPRB3
C
C  A second order boundary value problem with a parameter:
C
C  Y''+LAMBDA*EXP(Y)=0
C  Y(0)=0 Y'(1)=0
C
C  We use a finite difference approximation, with 21 equally
C  spaced nodes.  The value of Y at node (I) is X(I), for I=1 to 21.
C  X(22) is the value of LAMBDA.
C
C  We expect a limit point in LAMBDA at roughly LAMBDA=0.878.
C
C
C  This problem is solved six times:
C
C    With full storage user jacobian,
C    With full storage forward difference jacobian,
C    With full storage central difference jacobian.
C    With band storage user jacobian,
C    With band storage forward difference jacobian,
C    With band storage central difference jacobian.
C
      INTEGER   LIW
      INTEGER   LRW
      INTEGER   NVAR
C
      PARAMETER (NVAR=22)
      PARAMETER (LRW=29+(NVAR+6)*NVAR)
      PARAMETER (LIW=NVAR+29)
C
      EXTERNAL  BANSLV
      EXTERNAL  DENSLV
      EXTERNAL  FPBAND
      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   JAC
      INTEGER   LOUNIT
      CHARACTER NAME*12
      INTEGER   NIT
      REAL      RWORK(LRW)
      REAL      XR(NVAR)
C
      LOUNIT=6
C
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'PCPRB3'
      WRITE(LOUNIT,*)'PITCON sample program.'
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'Two point BVP, seeking limit point in lambda.'
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'This problem will be run SIX times.'
      WRITE(LOUNIT,*)' '
C
      ITRY=0
10    CONTINUE
C
      ITRY=ITRY+1
      WRITE(LOUNIT,*)' '
      WRITE(LOUNIT,*)'This is run number ',ITRY
      WRITE(LOUNIT,*)' '
C
      DO 20 I=1,LIW
        IWORK(I)=0
20      CONTINUE
      DO 30 I=1,LRW
       RWORK(I)=0.0
30     CONTINUE
      IERROR=0
C
      NIT=0
C
C  For runs 1 and 4, supply Jacobian.
C  For runs 2 and 5, use forward difference approximation,
C  For runs 3 and 6, use central difference approximation.
C
      IF(MOD(ITRY,3).EQ.1)JAC=0
      IF(MOD(ITRY,3).EQ.2)JAC=1
      IF(MOD(ITRY,3).EQ.0)JAC=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)=1    ; Cut down evaluations of jacobian on newton steps.
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)=*    ; Jacobian choice.  0=user, 1=forward, 2=central.
C
      IWORK(1)=0
      IWORK(2)=NVAR
      IWORK(3)=0
      IWORK(4)=1
      IWORK(5)=NVAR
      IWORK(6)=NVAR
      IWORK(7)=1
      IWORK(8)=LOUNIT
      IWORK(9)=JAC
C
C  Pass the lower and upper bandwidths of the Jacobian
C  in the integer parameter array IPAR.
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 80 I=1,NVAR
        XR(I)=0.0
80      CONTINUE
C
      IF(ITRY.LE.3)THEN
        WRITE(LOUNIT,*)'Use full solver DENSLV'
      ELSE
        WRITE(LOUNIT,*)'Use banded solver BANSLV'
        ENDIF
      IF(JAC.EQ.0)THEN
        WRITE(LOUNIT,*)'User supplies jacobian via subroutine.'
      ELSEIF(JAC.EQ.1)THEN
        WRITE(LOUNIT,*)'PITCON computes forward difference jacobian.'
      ELSEIF(JAC.EQ.2)THEN
        WRITE(LOUNIT,*)'PITCON computes center difference jacobian.'
        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
C  Take steps along the curve
C
      DO 90 I=1,50
        IF(ITRY.LE.3)THEN
          CALL PITCON(FPFULL,FPAR,FXBEND,IERROR,IPAR,IWORK,LIW,
     *    NVAR,RWORK,LRW,XR,DENSLV)
        ELSE
          CALL PITCON(FPBAND,FPAR,FXBEND,IERROR,IPAR,IWORK,LIW,
     *    NVAR,RWORK,LRW,XR,BANSLV)
          ENDIF
C
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,G14.6)')I,NAME,XR(NVAR)
C
        IF(IWORK(1).EQ.3)THEN
          WRITE(LOUNIT,*)' '
          NIT=NIT+1
          WRITE(LOUNIT,*)'Value of target point ',NIT
          WRITE(LOUNIT,'(1X,5G14.6)')(XR(J),J=1,NVAR)
          WRITE(LOUNIT,*)' '
          IF(NIT.GE.2)GO TO 100
          ENDIF
C
        IF(IERROR.NE.0)THEN
          WRITE(LOUNIT,*)'An error occurred, IERROR=',IERROR
          WRITE(LOUNIT,*)'We terminate this portion of the computation.'
          GO TO 100
          ENDIF
C
90      CONTINUE
C
100   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.6)GO TO 10
      STOP
      END
      SUBROUTINE FXBEND(NVAR,FPAR,IPAR,X,FX,IERROR)
C
C***********************************************************************
C
C  Nonlinear function evaluation for two point BVP.
C
C  Equation 1:  (boundary condition)
C
C  Y(0)=0  becomes FX(1)=X(1)
C
C  Equations 2 through NVAR-2 are the differential equation:
C
C  Y'' + LAMBDA*EXP(Y) = 0  becomes
C
C  FX(I)= X(I-1) -2*X(I) + X(I+1)
C        +X(NVAR)*EXP(X(I))*HSQ = 0
C
C  where HSQ is the square of the distance between nodes.
C
C  Equation NVAR-1:  (boundary condition)
C
C  Y'(1)=0  becomes  FX(NVAR-1)=X(NVAR-1)-X(NVAR-2)=0
C
      INTRINSIC EXP
C
      INTEGER   NVAR
C
      REAL      EXP
      REAL      FPAR(*)
      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
C
      INTRINSIC EXP
C
      INTEGER   NVAR
C
      REAL      EXP
      REAL      FPAR(*)
      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
      SUBROUTINE FPBAND(NVAR,FPAR,IPAR,X,FPRIME,IERROR)
C
C***********************************************************************
C
C  Band storage jacobian
C
      INTRINSIC EXP
C
      INTEGER   NVAR
C
      REAL      EXP
      REAL      FPAR(*)
      REAL      FPRIME(*)
      REAL      HSQ
      INTEGER   I
      INTEGER   IERROR
      INTEGER   INDX
      INTEGER   IPAR(2)
      INTEGER   ML
      INTEGER   MU
      INTEGER   NBAND
      REAL      X(NVAR)
C
      HSQ=NVAR-2
      HSQ=1.0/HSQ
      HSQ=HSQ*HSQ
      ML=IPAR(1)
      MU=IPAR(2)
      NBAND=2*ML+MU+1
      INDX=1+1*(NBAND-1)-ML
      FPRIME(INDX)=1.0
      DO 10 I=2,NVAR-2
        INDX=I+(I-1)*(NBAND-1)-ML
        FPRIME(INDX)=1.0
        INDX=I+I*(NBAND-1)-ML
        FPRIME(INDX)=-2.0+X(NVAR)*EXP(X(I))*HSQ
        INDX=I+(I+1)*(NBAND-1)-ML
        FPRIME(INDX)=1.0
        INDX=(NVAR-1)*NBAND+I
        FPRIME(INDX)=EXP(X(I))*HSQ
10      CONTINUE
      INDX=NVAR-1+(NVAR-2)*(NBAND-1)-ML
      FPRIME(INDX)=-1.0
      INDX=NVAR-1+(NVAR-1)*(NBAND-1)-ML
      FPRIME(INDX)=1.0
      RETURN
      END