      PARAMETER (MM = 5, NN = 3)
      DOUBLE PRECISION A(MM,NN), B(MM), X(NN), QRAUX(NN), WORK(NN), TOL
      INTEGER   JPVT(NN)
      DATA      B / 1.0D0, 2.3D0, 4.6D0, 3.1D0, 1.2D0/
C
C SET UP LEAST-SQUARES PROBLEM
C     QUADRATIC MODEL, EQUALLY-SPACED POINTS
C
      M = 5
      N = 3
      DO 10 I = 1,M
         A(I,1) = 1.0D0
         DO 10 J = 2,N
            A(I,J) = A(I,J-1)*I
10    CONTINUE
      TOL = 1.D-15
      ITASK = 1
      WRITE (*,*)   ' COEFFICIENT MATRIX'
      WRITE (*,800) ((A(I,J), J = 1,N), I = 1,M)
      WRITE (*,*)   ' RIGHT-HAND SIDE'
      WRITE (*,810) (B(I), I = 1,M)
C
C
C SOLVE LEAST-SQUARES PROBLEM
C
      CALL DQRLS (A, MM, M, N, TOL, KR, B, X, B, JPVT, QRAUX, WORK,
     *            ITASK, IND)
      IF (IND .NE. 0) THEN
         WRITE (*,*) ' ERROR CODE =', IND
         STOP
      END IF
      WRITE (*,*)    ' RANK OF MATRIX =', KR
C
C PRINT RESULTS
C
      WRITE (*,*)   ' PARAMETERS'
      WRITE (*,800) (X(J), J = 1,N)
      WRITE (*,*)   ' RESIDUALS'
      WRITE (*,810) (B(I), I = 1,M)
      WRITE (*,*)
      WRITE (*,*) 'REFERENCE RESULTS FROM IBM PC/AT'
      WRITE (*,*) ' RANK OF MATRIX =           3'
      WRITE (*,*) ' PARAMETERS'
      WRITE (*,*) '   -0.30200000E+01   0.44914286E+01  -0.72857143E+00'
      WRITE (*,*) ' RESIDUALS'
      WRITE (*,*) 
     *'   0.2571429E+00 -0.7485714E+00  0.7028571E+00 -0.1885714E+00
     *-0.2285714E-01'
C
      STOP
800   FORMAT (2X,3D17.8)
810   FORMAT (2X,5D15.7)
      END
