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