PROGRAM f08kgfe

!      F08KGF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : dgebrd, dgelqf, dgeqrf, dorglq, dorgqr, dormbr, &
                               f06qff, f06qhf, nag_wp, x04caf
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       REAL (KIND=nag_wp), PARAMETER   :: zero = 0.0E0_nag_wp
       INTEGER, PARAMETER              :: nin = 5, nout = 6
!      .. Local Scalars ..
       INTEGER                         :: i, ic, ifail, info, lda, ldpt, ldu,  &
                                          lwork, m, n
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: a(:,:), d(:), e(:), pt(:,:), tau(:), &
                                          taup(:), tauq(:), u(:,:), work(:)
!      .. Executable Statements ..
       WRITE (nout,*) 'F08KGF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
       DO ic = 1, 2
          READ (nin,*) m, n
          lda = m
          ldpt = n
          ldu = m
          lwork = 64*(m+n)
          ALLOCATE (a(lda,n),d(n),e(n-1),pt(ldpt,n),tau(n),taup(n),tauq(n), &
             u(ldu,n),work(lwork))

!         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 dgeqrf is f08aef
             CALL dgeqrf(m,n,a,lda,tau,work,lwork,info)

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

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

!            Copy R to PT (used as workspace)
             CALL f06qff('Upper',n,n,a,lda,pt,ldpt)

!            Set the strictly lower triangular part of R to zero
             CALL f06qhf('Lower',n-1,n-1,zero,zero,pt(2,1),ldpt)

!            Bidiagonalize R
!            The NAG name equivalent of dgebrd is f08kef
             CALL dgebrd(n,n,pt,ldpt,d,e,tauq,taup,work,lwork,info)

!            Update Q, storing the result in U
!            The NAG name equivalent of dormbr is f08kgf
             CALL dormbr('Q','Right','No transpose',m,n,n,pt,ldpt,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 x04caf('General',' ',m,n,u,ldu,'Example 1: matrix Q',ifail)

          ELSE

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

!            Copy A to PT
             CALL f06qff('Upper',m,n,a,lda,pt,ldpt)

!            Form Q explicitly, storing the result in PT
!            The NAG name equivalent of dorglq is f08ajf
             CALL dorglq(n,n,m,pt,ldpt,tau,work,lwork,info)

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

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

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

!            Update P**T, storing the result in PT
!            The NAG name equivalent of dormbr is f08kgf
             CALL dormbr('P','Left','Transpose',m,n,m,u,ldu,taup,pt,ldpt,work, &
                lwork,info)

!            Print bidiagonal form and matrix P**T

             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 x04caf('General',' ',m,n,pt,ldpt,'Example 2: matrix P**T', &
                ifail)

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

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