Program f08kgfe

!     F08KGF Example Program Text

!     Mark 28.3 Release. NAG Copyright 2022.

!     .. 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)

!         Make sign changes to standardise on positive diagonal elements.
          Do i = 1, n
            If (d(i)<zero) Then
              d(i) = -d(i)
              If (i<n) Then
                e(i) = -e(i)
              End If
              u(1:m,i) = -u(1:m,i)
            End If
          End Do

!         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,*) 'Superdiagonal'
          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)

!         Make sign changes to standardise on positive diagonal elements.
          Do i = 1, m
            If (d(i)<zero) Then
              d(i) = -d(i)
              If (i<m) Then
                e(i) = -e(i)
              End If
            End If
          End Do

!         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,*) 'Superdiagonal'
          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