*     F08ZEF Example Program Text
*     Mark 21 Release. NAG Copyright 2004.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          MMAX, NB, NMAX, PMAX
      PARAMETER        (MMAX=10,NB=64,NMAX=10,PMAX=10)
      INTEGER          LDA, LDB, LWORK
      PARAMETER        (LDA=NMAX,LDB=NMAX,LWORK=NB*(MMAX+PMAX))
      DOUBLE PRECISION ONE, ZERO
      PARAMETER        (ONE=1.0D0,ZERO=0.0D0)
*     .. Local Scalars ..
      DOUBLE PRECISION RNORM
      INTEGER          I, INFO, J, M, N, P
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,MMAX), B(LDB,PMAX), D(NMAX), TAUA(MMAX),
     +                 TAUB(MMAX+PMAX), WORK(LWORK), Y(PMAX)
*     .. External Functions ..
      DOUBLE PRECISION DNRM2
      EXTERNAL         DNRM2
*     .. External Subroutines ..
      EXTERNAL         DCOPY, DGEMV, DGGQRF, DORMQR, DORMRQ, DTRTRS,
     +                 F06FBF
*     .. Intrinsic Functions ..
      INTRINSIC        MAX, MIN
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F08ZEF Example Program Results'
      WRITE (NOUT,*)
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*) N, M, P
      IF (M.LE.MMAX .AND. N.LE.NMAX .AND. P.LE.PMAX .AND. M.LE.N .AND.
     +    N.LE.(M+P)) THEN
*
*        Read A, B and D from data file
*
         READ (NIN,*) ((A(I,J),J=1,M),I=1,N)
         READ (NIN,*) ((B(I,J),J=1,P),I=1,N)
         READ (NIN,*) (D(I),I=1,N)
*
*        Compute the generalized QR factorization of (A,B) as
*        A = Q*(R),   B = Q*(T11 T12)*Z
*              (0)          ( 0  T22)
*
         CALL DGGQRF(N,M,P,A,LDA,TAUA,B,LDB,TAUB,WORK,LWORK,INFO)
*
*        Compute c = (c1) = (Q**T)*d, storing the result in D
*                    (c2)
*
         CALL DORMQR('Left','Transpose',N,1,M,A,LDA,TAUA,D,N,WORK,LWORK,
     +               INFO)
*
*        Putting Z*y = w = (w1), set w1 = 0, storing the result in Y1
*                          (w2)
*
         CALL F06FBF(M+P-N,ZERO,Y,1)
*
         IF (N.GT.M) THEN
*
*           Copy c2 into Y2
*
            CALL DCOPY(N-M,D(M+1),1,Y(M+P-N+1),1)
*
*           Solve T22*w2 = c2 for w2, storing result in Y2
*
            CALL DTRTRS('Upper','No transpose','Non-unit',N-M,1,
     +                  B(M+1,M+P-N+1),LDB,Y(M+P-N+1),N-M,INFO)
*
            IF (INFO.GT.0) THEN
               WRITE (NOUT,*)
     +           'The upper triangular factor, T22, of B is singular, '
               WRITE (NOUT,*)
     +           'the least squares solution could not be computed'
               GO TO 20
            END IF
*
*           Compute estimate of the square root of the residual sum of
*           squares norm(y) = norm(w2)
*
            RNORM = DNRM2(N-M,Y(M+P-N+1),1)
*
*           Form c1 - T12*w2 in D
*
            CALL DGEMV('No transpose',M,N-M,-ONE,B(1,M+P-N+1),LDB,
     +                 Y(M+P-N+1),1,ONE,D,1)
         END IF
*
*        Solve R*x = c1 - T12*w2 for x
*
         CALL DTRTRS('Upper','No transpose','Non-unit',M,1,A,LDA,D,M,
     +               INFO)
*
         IF (INFO.GT.0) THEN
            WRITE (NOUT,*)
     +        'The upper triangular factor, R, of A is singular, '
            WRITE (NOUT,*)
     +        'the least squares solution could not be computed'
            GO TO 20
         END IF
*
*        Compute y = (Z**T)*w
*
         CALL DORMRQ('Left','Transpose',P,1,MIN(N,P),B(MAX(1,N-P+1),1),
     +               LDB,TAUB,Y,P,WORK,LWORK,INFO)
*
*        Print least squares solution x
*
         WRITE (NOUT,*) 'Generalized least squares solution'
         WRITE (NOUT,99999) (D(I),I=1,M)
*
*        Print residual vector y
*
         WRITE (NOUT,*)
         WRITE (NOUT,*) 'Residual vector'
         WRITE (NOUT,99998) (Y(I),I=1,P)
*
*        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 MMAX, NMAX or PMAX is too small, ',
     +     'and/or M.LT.N or N.GT.(M+P)'
      END IF
   20 CONTINUE
*
99999 FORMAT (1X,7F11.4)
99998 FORMAT (3X,1P,7E11.2)
      END