* E04GBF Example Program Text * Mark 15 Revised. NAG Copyright 1991. * .. Parameters .. INTEGER M, N, NT, LIW, LDFJAC, LDV, LW PARAMETER (M=15,N=3,NT=3,LIW=1,LDFJAC=M,LDV=N, + LW=7*N+M*N+2*M+N*N) INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) * .. Arrays in Common .. DOUBLE PRECISION T(M,NT), Y(M) * .. Local Scalars .. DOUBLE PRECISION ETA, FSUMSQ, STEPMX, XTOL INTEGER I, IFAIL, IPRINT, J, MAXCAL, NF, NITER * .. Local Arrays .. DOUBLE PRECISION FJAC(LDFJAC,N), FVEC(M), G(N), S(N), V(LDV,N), + W(LW), X(N) INTEGER IW(LIW) * .. External Functions .. DOUBLE PRECISION X02AJF EXTERNAL X02AJF * .. External Subroutines .. EXTERNAL E04GBF, E04HEV, E04YAF, LSQFUN, LSQGRD, LSQMON * .. Intrinsic Functions .. INTRINSIC SQRT * .. Common blocks .. COMMON Y, T * .. Executable Statements .. WRITE (NOUT,*) 'E04GBF Example Program Results' * Skip heading in data file READ (NIN,*) * Observations of TJ (J = 1, 2, 3) are held in T(I, J) * (I = 1, 2, . . . , 15) DO 20 I = 1, M READ (NIN,*) Y(I), (T(I,J),J=1,NT) 20 CONTINUE * Check LSQFUN by calling E04YAF at an arbitrary point X(1) = 0.19D0 X(2) = -1.34D0 X(3) = 0.88D0 IFAIL = 1 * CALL E04YAF(M,N,LSQFUN,X,FVEC,FJAC,LDFJAC,IW,LIW,W,LW,IFAIL) * IF (IFAIL.NE.0) THEN WRITE (NOUT,*) WRITE (NOUT,99999) ' ** E04YAF returned with IFAIL = ', IFAIL GO TO 40 END IF * Continue setting parameters for E04GBF * * Set IPRINT to 1 to obtain output from LSQMON at each iteration * IPRINT = -1 MAXCAL = 50*N * Since E04HEV is being used as LSQLIN, we set ETA to 0.9 ETA = 0.9D0 XTOL = 10.0D0*SQRT(X02AJF()) * We estimate that the minimum will be within 10 units of the * starting point STEPMX = 10.0D0 * Set up the starting point X(1) = 0.5D0 X(2) = 1.0D0 X(3) = 1.5D0 IFAIL = 1 * CALL E04GBF(M,N,E04HEV,LSQFUN,LSQMON,IPRINT,MAXCAL,ETA,XTOL, + STEPMX,X,FSUMSQ,FVEC,FJAC,LDFJAC,S,V,LDV,NITER,NF,IW, + LIW,W,LW,IFAIL) * IF (IFAIL.NE.0) THEN WRITE (NOUT,*) WRITE (NOUT,99999) ' ** E04GBF returned with IFAIL = ', IFAIL IF (IFAIL.LT.0 .OR. IFAIL.EQ.1) GO TO 40 END IF * WRITE (NOUT,*) WRITE (NOUT,99998) 'On exit, the sum of squares is', FSUMSQ WRITE (NOUT,99998) 'at the point', (X(J),J=1,N) CALL LSQGRD(M,N,FVEC,FJAC,LDFJAC,G) WRITE (NOUT,99997) 'The corresponding gradient is', (G(J),J=1,N) WRITE (NOUT,*) ' (machine dependent)' WRITE (NOUT,*) 'and the residuals are' WRITE (NOUT,99996) (FVEC(I),I=1,M) * 40 CONTINUE * 99999 FORMAT (1X,A,I5) 99998 FORMAT (1X,A,3F12.4) 99997 FORMAT (1X,A,1P,3E12.3) 99996 FORMAT (1X,1P,E9.1) END * SUBROUTINE LSQFUN(IFLAG,M,N,XC,FVEC,FJAC,LDFJAC,IW,LIW,W,LW) * Routine to evaluate the residuals and their 1st derivatives. * This routine is also suitable for use when E04FCV is used as * LSQLIN, since it can deal with IFLAG = 0 as well as IFLAG = 2. * .. Parameters .. INTEGER MDEC, NT PARAMETER (MDEC=15,NT=3) * .. Scalar Arguments .. INTEGER IFLAG, LDFJAC, LIW, LW, M, N * .. Array Arguments .. DOUBLE PRECISION FJAC(LDFJAC,N), FVEC(M), W(LW), XC(N) INTEGER IW(LIW) * .. Arrays in Common .. DOUBLE PRECISION T(MDEC,NT), Y(MDEC) * .. Local Scalars .. DOUBLE PRECISION DENOM, DUMMY INTEGER I * .. Common blocks .. COMMON Y, T * .. Executable Statements .. DO 20 I = 1, M DENOM = XC(2)*T(I,2) + XC(3)*T(I,3) FVEC(I) = XC(1) + T(I,1)/DENOM - Y(I) IF (IFLAG.EQ.0) GO TO 20 FJAC(I,1) = 1.0D0 DUMMY = -1.0D0/(DENOM*DENOM) FJAC(I,2) = T(I,1)*T(I,2)*DUMMY FJAC(I,3) = T(I,1)*T(I,3)*DUMMY 20 CONTINUE RETURN END * SUBROUTINE LSQMON(M,N,XC,FVEC,FJAC,LDFJAC,S,IGRADE,NITER,NF,IW, + LIW,W,LW) * Monitoring routine * .. Parameters .. INTEGER NDEC PARAMETER (NDEC=3) INTEGER NOUT PARAMETER (NOUT=6) * .. Scalar Arguments .. INTEGER IGRADE, LDFJAC, LIW, LW, M, N, NF, NITER * .. Array Arguments .. DOUBLE PRECISION FJAC(LDFJAC,N), FVEC(M), S(N), W(LW), XC(N) INTEGER IW(LIW) * .. Local Scalars .. DOUBLE PRECISION FSUMSQ, GTG INTEGER J * .. Local Arrays .. DOUBLE PRECISION G(NDEC) * .. External Functions .. DOUBLE PRECISION F06EAF EXTERNAL F06EAF * .. External Subroutines .. EXTERNAL LSQGRD * .. Executable Statements .. FSUMSQ = F06EAF(M,FVEC,1,FVEC,1) CALL LSQGRD(M,N,FVEC,FJAC,LDFJAC,G) GTG = F06EAF(N,G,1,G,1) WRITE (NOUT,*) WRITE (NOUT,*) + ' Itn F evals SUMSQ GTG Grade' WRITE (NOUT,99999) NITER, NF, FSUMSQ, GTG, IGRADE WRITE (NOUT,*) WRITE (NOUT,*) + ' X G Singular values' DO 20 J = 1, N WRITE (NOUT,99998) XC(J), G(J), S(J) 20 CONTINUE RETURN * 99999 FORMAT (1X,I4,6X,I5,6X,1P,E13.5,6X,1P,E9.1,6X,I3) 99998 FORMAT (1X,1P,E13.5,10X,1P,E9.1,10X,1P,E9.1) END * SUBROUTINE LSQGRD(M,N,FVEC,FJAC,LDFJAC,G) * Routine to evaluate gradient of the sum of squares * .. Scalar Arguments .. INTEGER LDFJAC, M, N * .. Array Arguments .. DOUBLE PRECISION FJAC(LDFJAC,N), FVEC(M), G(N) * .. Local Scalars .. DOUBLE PRECISION SUM INTEGER I, J * .. Executable Statements .. DO 40 J = 1, N SUM = 0.0D0 DO 20 I = 1, M SUM = SUM + FJAC(I,J)*FVEC(I) 20 CONTINUE G(J) = SUM + SUM 40 CONTINUE RETURN END