*     F08ZFF Example Program Text
*     Mark 21 Release. NAG Copyright 2004.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          PMAX, NB, NMAX, MMAX
      PARAMETER        (PMAX=10,NB=64,NMAX=10,MMAX=10)
      INTEGER          LDB, LDA, LWORK
      PARAMETER        (LDB=PMAX,LDA=MMAX,LWORK=NB*(PMAX+NMAX))
      DOUBLE PRECISION ONE
      PARAMETER        (ONE=1.0D0)
*     .. Local Scalars ..
      DOUBLE PRECISION RNORM
      INTEGER          I, INFO, J, M, N, P
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), B(LDB,NMAX), C(PMAX), D(MMAX),
     +                 TAUA(MIN(MMAX,NMAX)), TAUB(MIN(PMAX,NMAX)),
     +                 WORK(LWORK), X(NMAX)
*     .. External Functions ..
      DOUBLE PRECISION DNRM2
      EXTERNAL         DNRM2
*     .. External Subroutines ..
      EXTERNAL         DCOPY, DGEMV, DGGRQF, DORMQR, DORMRQ, DTRMV,
     +                 DTRTRS
*     .. Intrinsic Functions ..
      INTRINSIC        MIN
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F08ZFF Example Program Results'
      WRITE (NOUT,*)
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*) P, N, M
      IF (P.LE.PMAX .AND. N.LE.NMAX .AND. M.LE.MMAX .AND. M.LE.N .AND.
     +    N.LE.(P+M)) THEN
*
*        Read B, A, C and D from data file
*
         READ (NIN,*) ((B(I,J),J=1,N),I=1,P)
         READ (NIN,*) ((A(I,J),J=1,N),I=1,M)
         READ (NIN,*) (C(I),I=1,P)
         READ (NIN,*) (D(I),I=1,M)
*
*        Compute the generalized RQ factorization of (B,A) as
*        A = (0 R12)*Q,   B = Z*(T11 T12 T13)*Q, where R12, T11 and T22
*                               ( 0  T22 T23)
*        are upper triangular
*
         CALL DGGRQF(M,P,N,A,LDA,TAUA,B,LDB,TAUB,WORK,LWORK,INFO)
*
*        Compute (f1) = (Z**T)*c, storing the result in C
*                (f2)
*
         CALL DORMQR('Left','Transpose',P,1,MIN(P,N),B,LDB,TAUB,C,P,
     +               WORK,LWORK,INFO)
*
*        Putting Q*x = (y1), solve R12*w = d for w, storing result in D
*                      (w )
*
         CALL DTRTRS('Upper','No transpose','Non-unit',M,1,A(1,N-M+1),
     +               LDA,D,M,INFO)
*
         IF (INFO.GT.0) THEN
            WRITE (NOUT,*)
     +        'The upper triangular factor, R12, of A is singular, '
            WRITE (NOUT,*)
     +        'the least squares solution could not be computed'
            GO TO 40
         END IF
*
*        Form f1 - T1*w, T1 = (T12 T13), in C
*
         CALL DGEMV('No transpose',N-M,M,-ONE,B(1,N-M+1),LDB,D,1,ONE,C,
     +              1)
*
*        Solve T11*y1 = f1 - T1*w for y1, storing result in C
*
         CALL DTRTRS('Upper','No transpose','Non-unit',N-M,1,B,LDB,C,
     +               N-M,INFO)
*
         IF (INFO.GT.0) THEN
            WRITE (NOUT,*)
     +        'The upper triangular factor, T11, of B is singular, '
            WRITE (NOUT,*)
     +        'the least squares solution could not be computed'
            GO TO 40
         END IF
*
*        Copy y into X (first y1, then w)
*
         CALL DCOPY(N-M,C,1,X,1)
         CALL DCOPY(M,D,1,X(N-M+1),1)
*
*        Compute x = (Q**T)*y
*
         CALL DORMRQ('Left','Transpose',N,1,M,A,LDA,TAUA,X,N,WORK,LWORK,
     +               INFO)
*
*        Putting w = (y2), form f2 - T22*y2 - T23*y3
*                    (y3)
*
*        T22*y2
*
         CALL DTRMV('Upper','No transpose','Non-unit',MIN(P,N)-N+M,
     +              B(N-M+1,N-M+1),LDB,D,1)
*
*        f2 - T22*y2
*
         DO 20 I = 1, MIN(P,N) - N + M
            C(N-M+I) = C(N-M+I) - D(I)
   20    CONTINUE
         IF (P.LT.N) THEN
*
*           f2 - T22*y2 - T23*y3
*
            CALL DGEMV('No transpose',P-N+M,N-P,-ONE,B(N-M+1,P+1),LDB,
     +                 D(P-N+M+1),1,ONE,C(N-M+1),1)
         END IF
*
*        Compute estimate of the square root of the residual sum of
*        squares norm(r) = norm(f2 - T22*y2 - T23*y3)
*
         RNORM = DNRM2(P-(N-M),C(N-M+1),1)
*
*        Print least squares solution x
*
         WRITE (NOUT,*) 'Constrained least squares solution'
         WRITE (NOUT,99999) (X(I),I=1,N)
*
*        Print estimate of the square root of the residual sum of
*        squares
*
         WRITE (NOUT,*)
         WRITE (NOUT,*) 'Square root of the residual sum of squares'
         WRITE (NOUT,99998) RNORM
      ELSE
         WRITE (NOUT,*)
     +     'One or more of PMAX, NMAX or MMAX is too small, ',
     +     'and/or N.LT.M or N.GT.(P+M)'
      END IF
   40 CONTINUE
*
99999 FORMAT (1X,7F11.4)
99998 FORMAT (3X,1P,E11.2)
      END