
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**     DCOV_QQQ        **
C*************************
C
      SUBROUTINE DCOV_QQQ (CALFUN, M, N, X, FVEC, R, LDR, INFO, 
     +                     WA1, WA2, WA3, WA4)
C
C Reduced version of SLATEC library routine DCOV with IOPT=1 option only
C allowed. CALFUN (as for Laue library routines SSQMIN & LSQMIN and
C LSQMQL) replaces the normal FCN routine passed to DCOV which has additional 
C parameters not required by this version. Internal IFLAG cannot be reset via 
C user called function so INFO will not be returned as a negative value. 
C Messages also removed.
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 CALFUN
      DOUBLE PRECISION ONE,SIGMA,TEMP,ZERO,DENORM
      LOGICAL SING
      SAVE ZERO, ONE
      DATA ZERO/0.D0/,ONE/1.D0/
C
C-------------------------------------------------------------------------------
C
      SING=.FALSE.
      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 CALFUN (M,N,FVEC,X)
      TEMP=DENORM(M,FVEC)
      SIGMA=ONE
      IF (M.NE.N) SIGMA=TEMP*TEMP/(M-N)
C
C     CALCULATE THE JACOBIAN
C     CODE APPROXIMATES THE JACOBIAN
100   CALL DFDJC3_QQQ(CALFUN,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)
C
C     CHECK IF R IS SINGULAR.
      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 (SING) INFO=2
      RETURN
      END
