    Program f08ztfe

!     F08ZTF Example Program Text

!     Mark 24 Release. NAG Copyright 2012.

!     .. Use Statements ..
      Use nag_library, Only: dznrm2, nag_wp, zgemv, zggrqf, ztrmv, ztrtrs,     &
                             zunmqr, zunmrq
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Complex (Kind=nag_wp), Parameter :: one = (1.0E0_nag_wp,0.0E0_nag_wp)
      Integer, Parameter               :: nb = 64, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: rnorm
      Integer                          :: i, info, lda, ldb, lwork, m, n, p
!     .. Local Arrays ..
      Complex (Kind=nag_wp), Allocatable :: a(:,:), b(:,:), c(:), d(:),        &
                                            taua(:), taub(:), work(:), x(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: min
!     .. Executable Statements ..
      Write (nout,*) 'F08ZTF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) p, n, m
      lda = m
      ldb = p
      lwork = nb*(p+n)
      Allocate (a(lda,n),b(ldb,n),c(p),d(m),taua(min(m,n)),taub(min(p, &
        n)),work(lwork),x(n))

!     Read B, A, C and D from data file

      Read (nin,*)(b(i,1:n),i=1,p)
      Read (nin,*)(a(i,1:n),i=1,m)
      Read (nin,*) c(1:p)
      Read (nin,*) d(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
!     The NAG name equivalent of zggrqf is f08ztf
      Call zggrqf(m,p,n,a,lda,taua,b,ldb,taub,work,lwork,info)

!     Compute (f1) = (Z**H)*c, storing the result in C
!              (f2)
!     The NAG name equivalent of zunmqr is f08auf
      Call zunmqr('Left','Conjugate 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 )
!     The NAG name equivalent of ztrtrs is f07tsf
      Call ztrtrs('Upper','No transpose','Non-unit',m,1,a(1,n-m+1),lda,d,m, &
        info)

      If (info>0) Then
        Write (nout,*) 'The upper triangular factor, R12, of A is singular, '
        Write (nout,*) 'the least squares solution could not be computed'
      Else

!       Form f1 - T1*w, T1 = (T12 T13), in C

!       The NAG name equivalent of zgemv is f06saf
        Call zgemv('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
!       The NAG name equivalent of ztrtrs is f07tsf
        Call ztrtrs('Upper','No transpose','Non-unit',n-m,1,b,ldb,c,n-m,info)

        If (info>0) Then
          Write (nout,*) &
            'The upper triangular factor, T11, of B is singular, '
          Write (nout,*) 'the least squares solution could not be computed'
        Else

!         Copy y into X (first y1, then w)
          x(1:n-m) = c(1:n-m)
          x(n-m+1:n) = d(1:m)

!         Compute x = (Q**H)*y
!         The NAG name equivalent of zunmrq is f08cxf
          Call zunmrq('Left','Conjugate transpose',n,1,m,a,lda,taua,x,n,work, &
            lwork,info)

!         Putting w = (y2), form f2 - T22*y2 - T23*y3
!                     (y3)
!         T22*y2
!         The NAG name equivalent of ztrmv is f06sff
          Call ztrmv('Upper','No transpose','Non-unit',min(p,n)-n+m, &
            b(n-m+1,n-m+1),ldb,d,1)

!         f2 - T22*y2
          Do i = 1, min(p,n) - n + m
            c(n-m+i) = c(n-m+i) - d(i)
          End Do

          If (p<n) Then

!           f2 - T22*y2 - T23*y3
!           The NAG name equivalent of zgemv is f06saf
            Call zgemv('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)
!         The NAG name equivalent of dznrm2 is f06jjf
          rnorm = dznrm2(p-(n-m),c(n-m+1),1)

!         Print least squares solution x
          Write (nout,*) 'Constrained least squares solution'
          Write (nout,99999) x(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
        End If
      End If

99999 Format (4(' (',F7.4,',',F7.4,')':))
99998 Format (1X,1P,E10.2)
    End Program f08ztfe
