relibc_openlibm/slatec/dcov.f

*DECK DCOV
      SUBROUTINE DCOV (FCN, IOPT, M, N, X, FVEC, R, LDR, INFO, WA1, WA2,
     +   WA3, WA4)
C***BEGIN PROLOGUE  DCOV
C***PURPOSE  Calculate the covariance matrix for a nonlinear data
C            fitting problem.  It is intended to be used after a
C            successful return from either DNLS1 or DNLS1E.
C***LIBRARY   SLATEC
C***CATEGORY  K1B1
C***TYPE      DOUBLE PRECISION (SCOV-S, DCOV-D)
C***KEYWORDS  COVARIANCE MATRIX, NONLINEAR DATA FITTING,
C             NONLINEAR LEAST SQUARES
C***AUTHOR  Hiebert, K. L., (SNLA)
C***DESCRIPTION
C
C  1. Purpose.
C
C     DCOV calculates the covariance matrix for a nonlinear data
C     fitting problem.  It is intended to be used after a
C     successful return from either DNLS1 or DNLS1E. DCOV
C     and DNLS1 (and DNLS1E) have compatible parameters.  The
C     required external subroutine, FCN, is the same
C     for all three codes, DCOV, DNLS1, and DNLS1E.
C
C  2. Subroutine and Type Statements.
C
C     SUBROUTINE DCOV(FCN,IOPT,M,N,X,FVEC,R,LDR,INFO,
C                     WA1,WA2,WA3,WA4)
C     INTEGER IOPT,M,N,LDR,INFO
C     DOUBLE PRECISION X(N),FVEC(M),R(LDR,N),WA1(N),WA2(N),WA3(N),WA4(M)
C     EXTERNAL FCN
C
C  3. Parameters. All TYPE REAL parameters are DOUBLE PRECISION
C
C      FCN is the name of the user-supplied subroutine which calculates
C         the functions.  If the user wants to supply the Jacobian
C         (IOPT=2 or 3), then FCN must be written to calculate the
C         Jacobian, as well as the functions.  See the explanation
C         of the IOPT argument below.
C         If the user wants the iterates printed in DNLS1 or DNLS1E,
C         then FCN must do the printing.  See the explanation of NPRINT
C         in DNLS1 or DNLS1E.  FCN must be declared in an EXTERNAL
C         statement in the calling program and should be written as
C         follows.
C
C         SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
C         INTEGER IFLAG,LDFJAC,M,N
C         DOUBLE PRECISION X(N),FVEC(M)
C         ----------
C         FJAC and LDFJAC may be ignored       , if IOPT=1.
C         DOUBLE PRECISION FJAC(LDFJAC,N)      , if IOPT=2.
C         DOUBLE PRECISION FJAC(N)             , if IOPT=3.
C         ----------
C           If IFLAG=0, the values in X and FVEC are available
C           for printing in DNLS1 or DNLS1E.
C           IFLAG will never be zero when FCN is called by DCOV.
C           The values of X and FVEC must not be changed.
C         RETURN
C         ----------
C           If IFLAG=1, calculate the functions at X and return
C           this vector in FVEC.
C         RETURN
C         ----------
C           If IFLAG=2, calculate the full Jacobian at X and return
C           this matrix in FJAC.  Note that IFLAG will never be 2 unless
C           IOPT=2.  FVEC contains the function values at X and must
C           not be altered.  FJAC(I,J) must be set to the derivative
C           of FVEC(I) with respect to X(J).
C         RETURN
C         ----------
C           If IFLAG=3, calculate the LDFJAC-th row of the Jacobian
C           and return this vector in FJAC.  Note that IFLAG will
C           never be 3 unless IOPT=3.  FJAC(J) must be set to
C           the derivative of FVEC(LDFJAC) with respect to X(J).
C         RETURN
C         ----------
C         END
C
C
C         The value of IFLAG should not be changed by FCN unless the
C         user wants to terminate execution of DCOV.  In this case, set
C         IFLAG to a negative integer.
C
C
C       IOPT is an input variable which specifies how the Jacobian will
C         be calculated.  If IOPT=2 or 3, then the user must supply the
C         Jacobian, as well as the function values, through the
C         subroutine FCN.  If IOPT=2, the user supplies the full
C         Jacobian with one call to FCN.  If IOPT=3, the user supplies
C         one row of the Jacobian with each call.  (In this manner,
C         storage can be saved because the full Jacobian is not stored.)
C         If IOPT=1, the code will approximate the Jacobian by forward
C         differencing.
C
C       M is a positive integer input variable set to the number of
C         functions.
C
C       N is a positive integer input variable set to the number of
C         variables.  N must not exceed M.
C
C       X is an array of length N.  On input X must contain the value
C         at which the covariance matrix is to be evaluated.  This is
C         usually the value for X returned from a successful run of
C         DNLS1 (or DNLS1E).  The value of X will not be changed.
C
C    FVEC is an output array of length M which contains the functions
C         evaluated at X.
C
C       R is an output array.  For IOPT=1 and 2, R is an M by N array.
C         For IOPT=3, R is an N by N array.  On output, if INFO=1,
C         the upper N by N submatrix of R contains the covariance
C         matrix evaluated at X.
C
C     LDR is a positive integer input variable which specifies
C         the leading dimension of the array R.  For IOPT=1 and 2,
C         LDR must not be less than M.  For IOPT=3, LDR must not
C         be less than N.
C
C    INFO is an integer output variable.  If the user has terminated
C         execution, INFO is set to the (negative) value of IFLAG.  See
C         description of FCN. Otherwise, INFO is set as follows.
C
C         INFO = 0 Improper input parameters (M.LE.0 or N.LE.0).
C
C         INFO = 1 Successful return.  The covariance matrix has been
C                  calculated and stored in the upper N by N
C                  submatrix of R.
C
C         INFO = 2 The Jacobian matrix is singular for the input value
C                  of X.  The covariance matrix cannot be calculated.
C                  The upper N by N submatrix of R contains the QR
C                  factorization of the Jacobian (probably not of
C                  interest to the user).
C
C WA1,WA2 are work arrays of length N.
C and WA3
C
C     WA4 is a work array of length M.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DENORM, DFDJC3, DQRFAC, DWUPDT, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   810522  DATE WRITTEN
C   890831  Modified array declarations.  (WRB)
C   891006  Cosmetic changes to prologue.  (WRB)
C   891006  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900510  Fixed an error message.  (RWC)
C***END PROLOGUE  DCOV
C
C     REVISED 850601-1100
C     REVISED YYMMDD HHMM
C
      INTEGER I,IDUM,IFLAG,INFO,IOPT,J,K,KP1,LDR,M,N,NM1,NROW
      DOUBLE PRECISION X(*),R(LDR,*),FVEC(*),WA1(*),WA2(*),WA3(*),
     1 WA4(*)
      EXTERNAL FCN
      DOUBLE PRECISION ONE,SIGMA,TEMP,ZERO,DENORM
      LOGICAL SING
      SAVE ZERO, ONE
      DATA ZERO/0.D0/,ONE/1.D0/
C***FIRST EXECUTABLE STATEMENT  DCOV
      SING=.FALSE.
      IFLAG=0
      IF (M.LE.0 .OR. N.LE.0) GO TO 300
C
C     CALCULATE SIGMA = (SUM OF THE SQUARED RESIDUALS) / (M-N)
      IFLAG=1
      CALL FCN(IFLAG,M,N,X,FVEC,R,LDR)
      IF (IFLAG.LT.0) GO TO 300
      TEMP=DENORM(M,FVEC)
      SIGMA=ONE
      IF (M.NE.N) SIGMA=TEMP*TEMP/(M-N)
C
C     CALCULATE THE JACOBIAN
      IF (IOPT.EQ.3) GO TO 200
C
C     STORE THE FULL JACOBIAN USING M*N STORAGE
      IF (IOPT.EQ.1) GO TO 100
C
C     USER SUPPLIES THE JACOBIAN
      IFLAG=2
      CALL FCN(IFLAG,M,N,X,FVEC,R,LDR)
      GO TO 110
C
C     CODE APPROXIMATES THE JACOBIAN
100   CALL DFDJC3(FCN,M,N,X,FVEC,R,LDR,IFLAG,ZERO,WA4)
110   IF (IFLAG.LT.0) GO TO 300
C
C     COMPUTE THE QR DECOMPOSITION
      CALL DQRFAC(M,N,R,LDR,.FALSE.,IDUM,1,WA1,WA1,WA1)
      DO 120 I=1,N
120   R(I,I)=WA1(I)
      GO TO 225
C
C     COMPUTE THE QR FACTORIZATION OF THE JACOBIAN MATRIX CALCULATED ONE
C     ROW AT A TIME AND STORED IN THE UPPER TRIANGLE OF R.
C     ( (Q TRANSPOSE)*FVEC IS ALSO CALCULATED BUT NOT USED.)
200   CONTINUE
      DO 210 J=1,N
      WA2(J)=ZERO
      DO 205 I=1,N
      R(I,J)=ZERO
205   CONTINUE
210   CONTINUE
      IFLAG=3
      DO 220 I=1,M
      NROW = I
      CALL FCN(IFLAG,M,N,X,FVEC,WA1,NROW)
      IF (IFLAG.LT.0) GO TO 300
      TEMP=FVEC(I)
      CALL DWUPDT(N,R,LDR,WA1,WA2,TEMP,WA3,WA4)
220   CONTINUE
C
C     CHECK IF R IS SINGULAR.
225   CONTINUE
      DO 230 I=1,N
      IF (R(I,I).EQ.ZERO) SING=.TRUE.
230   CONTINUE
      IF (SING) GO TO 300
C
C     R IS UPPER TRIANGULAR.  CALCULATE (R TRANSPOSE) INVERSE AND STORE
C     IN THE UPPER TRIANGLE OF R.
      IF (N.EQ.1) GO TO 275
      NM1=N-1
      DO 270 K=1,NM1
C
C     INITIALIZE THE RIGHT-HAND SIDE (WA1(*)) AS THE K-TH COLUMN OF THE
C     IDENTITY MATRIX.
      DO 240 I=1,N
      WA1(I)=ZERO
240   CONTINUE
      WA1(K)=ONE
C
      R(K,K)=WA1(K)/R(K,K)
      KP1=K+1
      DO 260 I=KP1,N
C
C     SUBTRACT R(K,I-1)*R(I-1,*) FROM THE RIGHT-HAND SIDE, WA1(*).
      DO 250 J=I,N
      WA1(J)=WA1(J)-R(K,I-1)*R(I-1,J)
250   CONTINUE
      R(K,I)=WA1(I)/R(I,I)
260   CONTINUE
270   CONTINUE
275   R(N,N)=ONE/R(N,N)
C
C     CALCULATE R-INVERSE * (R TRANSPOSE) INVERSE AND STORE IN THE UPPER
C     TRIANGLE OF R.
      DO 290 I=1,N
      DO 290 J=I,N
      TEMP=ZERO
      DO 280 K=J,N
      TEMP=TEMP+R(I,K)*R(J,K)
280   CONTINUE
      R(I,J)=TEMP*SIGMA
290   CONTINUE
      INFO=1
C
300   CONTINUE
      IF (M.LE.0 .OR. N.LE.0) INFO=0
      IF (IFLAG.LT.0) INFO=IFLAG
      IF (SING) INFO=2
      IF (INFO .LT. 0) CALL XERMSG ('SLATEC', 'DCOV',
     +   'EXECUTION TERMINATED BECAUSE USER SET IFLAG NEGATIVE.', 1, 1)
      IF (INFO .EQ. 0) CALL XERMSG ('SLATEC', 'DCOV',
     +   'INVALID INPUT PARAMETER.', 2, 1)
      IF (INFO .EQ. 2) CALL XERMSG ('SLATEC', 'DCOV',
     +   'SINGULAR JACOBIAN MATRIX, COVARIANCE MATRIX CANNOT BE ' //
     +   'CALCULATED.', 1, 1)
      RETURN
      END