PROGRAM f08kufe

!      F08KUF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : f06tff, f06thf, nag_wp, x04dbf, zgebrd, zgelqf, &
                               zgeqrf, zunglq, zungqr, zunmbr
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       COMPLEX (KIND=nag_wp), PARAMETER :: zero = (0.0E0_nag_wp,0.0E0_nag_wp)
       INTEGER, PARAMETER              :: nin = 5, nout = 6
!      .. Local Scalars ..
       INTEGER                         :: i, ic, ifail, info, lda, ldph, ldu,  &
                                          lwork, m, n
!      .. Local Arrays ..
       COMPLEX (KIND=nag_wp), ALLOCATABLE :: a(:,:), ph(:,:), tau(:), taup(:), &
                                             tauq(:), u(:,:), work(:)
       REAL (KIND=nag_wp), ALLOCATABLE :: d(:), e(:)
       CHARACTER (1)                   :: clabs(1), rlabs(1)
!      .. Executable Statements ..
       WRITE (nout,*) 'F08KUF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
       DO ic = 1, 2
          READ (nin,*) m, n
          lda = m
          ldph = n
          ldu = m
          lwork = 64*(m+n)
          ALLOCATE (a(lda,n),ph(ldph,n),tau(n),taup(n),tauq(n),u(ldu,n), &
             work(lwork),d(n),e(n-1))

!         Read A from data file

          READ (nin,*) (a(i,1:n),i=1,m)

          IF (m>=n) THEN

!            Compute the QR factorization of A
!            The NAG name equivalent of zgeqrf is f08asf
             CALL zgeqrf(m,n,a,lda,tau,work,lwork,info)

!            Copy A to U
             CALL f06tff('Lower',m,n,a,lda,u,ldu)

!            Form Q explicitly, storing the result in U
!            The NAG name equivalent of zungqr is f08atf
             CALL zungqr(m,n,n,u,ldu,tau,work,lwork,info)

!            Copy R to PH (used as workspace)
             CALL f06tff('Upper',n,n,a,lda,ph,ldph)

!            Set the strictly lower triangular part of R to zero
             CALL f06thf('Lower',n-1,n-1,zero,zero,ph(2,1),ldph)

!            Bidiagonalize R
!            The NAG name equivalent of zgebrd is f08ksf
             CALL zgebrd(n,n,ph,ldph,d,e,tauq,taup,work,lwork,info)

!            Update Q, storing the result in U
!            The NAG name equivalent of zunmbr is f08kuf
             CALL zunmbr('Q','Right','No transpose',m,n,n,ph,ldph,tauq,u,ldu, &
                work,lwork,info)

!            Print bidiagonal form and matrix Q

             WRITE (nout,*)
             WRITE (nout,*) 'Example 1: bidiagonal matrix B'
             WRITE (nout,*) 'Diagonal'
             WRITE (nout,99999) d(1:n)
             WRITE (nout,*) 'Super-diagonal'
             WRITE (nout,99999) e(1:n-1)
             WRITE (nout,*)
             FLUSH (nout)

!            ifail: behaviour on error exit
!                   =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
             ifail = 0
             CALL x04dbf('General',' ',m,n,u,ldu,'Bracketed','F7.4', &
                'Example 1: matrix Q','Integer',rlabs,'Integer',clabs,80,0, &
                ifail)

          ELSE

!            Compute the LQ factorization of A
!            The NAG name equivalent of zgelqf is f08avf
             CALL zgelqf(m,n,a,lda,tau,work,lwork,info)

!            Copy A to PH
             CALL f06tff('Upper',m,n,a,lda,ph,ldph)

!            Form Q explicitly, storing the result in PH
!            The NAG name equivalent of zunglq is f08awf
             CALL zunglq(n,n,m,ph,ldph,tau,work,lwork,info)

!            Copy L to U (used as workspace)
             CALL f06tff('Lower',m,m,a,lda,u,ldu)

!            Set the strictly upper triangular part of L to zero
             CALL f06thf('Upper',m-1,m-1,zero,zero,u(1,2),ldu)

!            Bidiagonalize L
!            The NAG name equivalent of zgebrd is f08ksf
             CALL zgebrd(m,m,u,ldu,d,e,tauq,taup,work,lwork,info)

!            Update P**H, storing the result in PH
!            The NAG name equivalent of zunmbr is f08kuf
             CALL zunmbr('P','Left','Conjugate transpose',m,n,m,u,ldu,taup,ph, &
                ldph,work,lwork,info)

!            Print bidiagonal form and matrix P**H

             WRITE (nout,*)
             WRITE (nout,*) 'Example 2: bidiagonal matrix B'
             WRITE (nout,*) 'Diagonal'
             WRITE (nout,99999) d(1:m)
             WRITE (nout,*) 'Super-diagonal'
             WRITE (nout,99999) e(1:m-1)
             WRITE (nout,*)
             FLUSH (nout)

             ifail = 0
             CALL x04dbf('General',' ',m,n,ph,ldph,'Bracketed','F7.4', &
                'Example 2: matrix P**H','Integer',rlabs,'Integer',clabs,80,0, &
                ifail)

          END IF
          DEALLOCATE (a,ph,tau,taup,tauq,u,work,d,e)
       END DO

99999  FORMAT (3X,(8F8.4))
    END PROGRAM f08kufe