*     F08ZTF 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=MMAX,LDB=PMAX,LWORK=NB*(MMAX+NMAX))
      COMPLEX *16      ONE
      PARAMETER        (ONE=(1.0D0,0.0D0))
*     .. Local Scalars ..
      DOUBLE PRECISION RNORM
      INTEGER          I, INFO, J, M, N, P
*     .. Local Arrays ..
      COMPLEX *16      A(LDA,NMAX), B(LDB,NMAX), C(MMAX), D(PMAX),
     +                 TAUA(MMAX+NMAX), TAUB(PMAX), WORK(LWORK), X(NMAX)
*     .. External Functions ..
      DOUBLE PRECISION DZNRM2
      EXTERNAL         DZNRM2
*     .. External Subroutines ..
      EXTERNAL         ZCOPY, ZGEMV, ZGGRQF, ZTRMV, ZTRTRS, ZUNMQR,
     +                 ZUNMRQ
*     .. Intrinsic Functions ..
      INTRINSIC        MIN
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F08ZTF Example Program Results'
      WRITE (NOUT,*)
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*) M, N, P
      IF (M.LE.MMAX .AND. N.LE.NMAX .AND. P.LE.PMAX .AND. P.LE.N .AND.
     +    N.LE.(M+P)) THEN
*
*        Read A, B, c and d from data file
*
         READ (NIN,*) ((A(I,J),J=1,N),I=1,M)
         READ (NIN,*) ((B(I,J),J=1,N),I=1,P)
         READ (NIN,*) (C(I),I=1,M)
         READ (NIN,*) (D(I),I=1,P)
*
*        Compute the generalized RQ factorization of (A,B) as
*        B = (0 R12)*Q,   A = Z*(T11 T12 T13)*Q, where R12, T11 and T22
*                               ( 0  T22 T23)
*        are upper triangular
*
         CALL ZGGRQF(P,M,N,B,LDB,TAUB,A,LDA,TAUA,WORK,LWORK,INFO)
*
*        Compute (f1) = (Z**H)*c, storing the result in C
*                (f2)
*
         CALL ZUNMQR('Left','Conjugate transpose',M,1,MIN(M,N),A,LDA,
     +               TAUA,C,M,WORK,LWORK,INFO)
*
*        Putting Q*x = (y1), solve R12*w = d for w, storing result in D
*                      (w )
*
         CALL ZTRTRS('Upper','No transpose','Non-unit',P,1,B(1,N-P+1),
     +               LDB,D,P,INFO)
*
         IF (INFO.GT.0) THEN
            WRITE (NOUT,*)
     +        'The upper triangular factor, R12, of B 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 ZGEMV('No transpose',N-P,P,-ONE,A(1,N-P+1),LDA,D,1,ONE,C,
     +              1)
*
*        Solve T11*y1 = f1 - T1*w for y1, storing result in C
*
         CALL ZTRTRS('Upper','No transpose','Non-unit',N-P,1,A,LDA,C,
     +               N-P,INFO)
*
         IF (INFO.GT.0) THEN
            WRITE (NOUT,*)
     +        'The upper triangular factor, T11, of A 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 ZCOPY(N-P,C,1,X,1)
         CALL ZCOPY(P,D,1,X(N-P+1),1)
*
*        Compute x = (Q**H)*y
*
         CALL ZUNMRQ('Left','Conjugate transpose',N,1,P,B,LDB,TAUB,X,N,
     +               WORK,LWORK,INFO)
*
*        Putting w = (y2), form f2 - T22*y2 - T23*y3
*                    (y3)
*
*        T22*y2
*
         CALL ZTRMV('Upper','No transpose','Non-unit',MIN(M,N)-N+P,
     +              A(N-P+1,N-P+1),LDA,D,1)
*
*        f2 - T22*y2
*
         DO 20 I = 1, MIN(M,N) - N + P
            C(N-P+I) = C(N-P+I) - D(I)
   20    CONTINUE
         IF (M.LT.N) THEN
*
*           f2 - T22*y2 - T23*y3
*
            CALL ZGEMV('No transpose',M-N+P,N-M,-ONE,A(N-P+1,M+1),LDA,
     +                 D(M-N+P+1),1,ONE,C(N-P+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 = DZNRM2(M-(N-P),C(N-P+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 MMAX, NMAX or PMAX is too small, ',
     +     'and/or N.LT.P or N.GT.(M+P)'
      END IF
   40 CONTINUE
      STOP
*
99999 FORMAT (4(' (',F7.4,',',F7.4,')',:))
99998 FORMAT (1X,1P,E10.2)
      END