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