* E04YBF Example Program Text * Mark 14 Revised. NAG Copyright 1989. * .. Parameters .. INTEGER MDEC, NDEC, LDFJAC, LB, LIW, LW PARAMETER (MDEC=15,NDEC=3,LDFJAC=MDEC,LB=NDEC*(NDEC+1)/2, + LIW=1,LW=5*NDEC+MDEC+MDEC*NDEC+NDEC*(NDEC-1)/2) INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) * .. Arrays in Common .. DOUBLE PRECISION T(MDEC,NDEC), Y(MDEC) * .. Local Scalars .. INTEGER I, IFAIL, J, K, M, N * .. Local Arrays .. DOUBLE PRECISION B(LB), FJAC(LDFJAC,NDEC), FVEC(MDEC), W(LW), + X(NDEC) INTEGER IW(LIW) * .. External Subroutines .. EXTERNAL E04YAF, E04YBF, LSQFUN, LSQHES * .. Common blocks .. COMMON Y, T * .. Executable Statements .. WRITE (NOUT,*) 'E04YBF Example Program Results' * Skip heading in data file READ (NIN,*) M = MDEC N = NDEC * 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,N) 20 CONTINUE * Set up an arbitrary point at which to check the derivatives X(1) = 0.19D0 X(2) = -1.34D0 X(3) = 0.88D0 * Check the 1st derivatives IFAIL = 1 * CALL E04YAF(M,N,LSQFUN,X,FVEC,FJAC,LDFJAC,IW,LIW,W,LW,IFAIL) IF (IFAIL.NE.0) THEN WRITE (NOUT,99998) IFAIL GO TO 60 END IF * WRITE (NOUT,*) WRITE (NOUT,*) 'The test point is' WRITE (NOUT,99999) (X(J),J=1,N) * * Check the evaluation of B IFAIL = 1 * CALL E04YBF(M,N,LSQFUN,LSQHES,X,FVEC,FJAC,LDFJAC,B,LB,IW,LIW,W,LW, + IFAIL) * WRITE (NOUT,*) IF (IFAIL.EQ.1) THEN WRITE (NOUT,*) 'A parameter is outside its expected range' ELSE IF (IFAIL.GE.0) THEN IF (IFAIL.EQ.0) THEN WRITE (NOUT,*) + 'The matrix B is consistent with 1st derivatives' ELSE IF (IFAIL.EQ.2) THEN WRITE (NOUT,*) + 'Probable error in calculation of the matrix B' END IF WRITE (NOUT,*) WRITE (NOUT,*) 'At the test point, LSQFUN gives' WRITE (NOUT,*) WRITE (NOUT,*) + ' Residuals 1st derivatives' WRITE (NOUT,99997) (FVEC(I),(FJAC(I,J),J=1,N),I=1,M) WRITE (NOUT,*) WRITE (NOUT,*) + 'and LSQHES gives the lower triangle of the matrix B' WRITE (NOUT,*) K = 1 DO 40 I = 1, N WRITE (NOUT,99997) (B(J),J=K,K+I-1) K = K + I 40 CONTINUE ELSE IF (IFAIL.EQ.-399) THEN WRITE (NOUT,99996) IFAIL ELSE WRITE (NOUT,99995) 'IFLAG was set to ', IFAIL, + 'in LSQFUN or LSQHES' END IF 60 CONTINUE * 99999 FORMAT (1X,4F10.5) 99998 FORMAT (1X,/1X,' ** E04YAF returned with IFAIL = ',I5) 99997 FORMAT (1X,1P,4E15.3) 99996 FORMAT (1X,/1X,' ** E04YBF returned with IFAIL = ',I5) 99995 FORMAT (1X,A,I3,A) END * SUBROUTINE LSQFUN(IFLAG,M,N,XC,FVEC,FJAC,LDFJAC,IW,LIW,W,LW) * Routine to evaluate the residuals and their 1st derivatives * .. Parameters .. INTEGER MDEC, NDEC PARAMETER (MDEC=15,NDEC=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,NDEC), 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) 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 LSQHES(IFLAG,M,N,FVEC,XC,B,LB,IW,LIW,W,LW) * Routine to compute the lower triangle of the matrix B * (stored by rows in the array B) * .. Parameters .. INTEGER MDEC, NDEC PARAMETER (MDEC=15,NDEC=3) * .. Scalar Arguments .. INTEGER IFLAG, LB, LIW, LW, M, N * .. Array Arguments .. DOUBLE PRECISION B(LB), FVEC(M), W(LW), XC(N) INTEGER IW(LIW) * .. Arrays in Common .. DOUBLE PRECISION T(MDEC,NDEC), Y(MDEC) * .. Local Scalars .. DOUBLE PRECISION DUMMY, SUM22, SUM32, SUM33 INTEGER I * .. Common blocks .. COMMON Y, T * .. Executable Statements .. B(1) = 0.0D0 B(2) = 0.0D0 SUM22 = 0.0D0 SUM32 = 0.0D0 SUM33 = 0.0D0 DO 20 I = 1, M DUMMY = 2.0D0*T(I,1)/(XC(2)*T(I,2)+XC(3)*T(I,3))**3 SUM22 = SUM22 + FVEC(I)*DUMMY*T(I,2)**2 SUM32 = SUM32 + FVEC(I)*DUMMY*T(I,2)*T(I,3) SUM33 = SUM33 + FVEC(I)*DUMMY*T(I,3)**2 20 CONTINUE B(3) = SUM22 B(4) = 0.0D0 B(5) = SUM32 B(6) = SUM33 RETURN END