Program f10dafe

!     F10DAF Example Program
!     Mark 28.3 Release. NAG Copyright 2022.

!     .. Use Statements ..
      Use nag_library, Only: dgemm, dgemv, dgeqrf, dgesvd, dnrm2, dorgqr,      &
                             f10daf, g05kff, g05skf, nag_wp, x01aaf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: zero = 0.0_nag_wp
      Integer, Parameter               :: lseed = 1, nb = 64, nin = 5,         &
                                          nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: err_bound, pi
      Integer                          :: genid, i, ifail, info, j, k, lda,    &
                                          ldu, ldvt, lstate, lwork, m, n, r,   &
                                          subid
      Character (1)                    :: trans
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), imqq(:,:), imqqa(:,:),       &
                                          omega(:), pa(:,:), q(:,:), res(:),   &
                                          res_norm(:), s(:), tau(:), work(:),  &
                                          y(:)
      Real (Kind=nag_wp)               :: dummy(1,1), u(1,1), vt(1,1)
      Integer                          :: seed(lseed)
      Integer, Allocatable             :: state(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max, maxval, nint, sqrt
!     .. Executable Statements ..
      Write (nout,*) 'F10DAF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) m, n
      lda = m
      Allocate (a(lda,n))

!     Read the m by n matrix A from data file
      Read (nin,*)(a(i,1:n),i=1,m)

!     Read in the base generator information and seed
      Read (nin,*) genid, subid, seed(1)

!     Read in trans, k - arguments to random projection routine
      Read (nin,*) trans, k

!     Read in r - parameter of error estimation
      Read (nin,*) r

      Write (nout,*) ' The matrix A'
      Write (nout,99999)(a(i,1:n),i=1,m)
99999 Format (6(3X,F5.2))

!     Initialise random number generator (state argument)

!     Initial call to initializer to get size of STATE array
      lstate = 0
      Allocate (state(lstate))
      ifail = 0
      Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!     Reallocate STATE
      Deallocate (state)
      Allocate (state(lstate))

!     Initialize the generator to a repeatable sequence
      ifail = 0
      Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!     Perform random projection
      Allocate (pa(m,n))
      ifail = 0
      pa(:,:) = a
      Call f10daf(trans,m,n,pa,lda,k,state,ifail)

!     Perform QR decomposition on random projection

!     QR decomposition of pa
      lwork = k*nb
      Allocate (tau(k),work(lwork),q(m,k))
      Call dgeqrf(m,k,pa,m,tau,work,lwork,info)
      q(:,:) = pa(1:m,1:k)
      Call dorgqr(m,k,k,q,m,tau,work,lwork,info)

!     Calculate true error as
!       ||(I - Q Q') A|| = largest singular value of (I - Q Q') A
      Allocate (imqq(m,m))

!     (I - Q Q')
      imqq(:,:) = zero
      Do i = 1, m
        imqq(i,i) = one
      End Do
      Call dgemm('N','T',m,m,k,-one,q,m,q,m,one,imqq,m)

!     (I - Q Q') A
      Allocate (imqqa(m,n))
      ldu = 1
      ldvt = 1
      Call dgemm('N','N',m,n,m,one,imqq,m,a,m,zero,imqqa,m)

      Allocate (s(m))
!     Use routine workspace query to get optimal workspace.
      lwork = -1
      Call dgesvd('N','N',m,n,imqqa,lda,s,u,ldu,vt,ldvt,dummy,lwork,info)

!     Make sure that there is enough workspace for block size nb.
      lwork = max(m+4*n+nb*(m+n),nint(dummy(1,1)))
      Deallocate (work)
      Allocate (work(lwork))

!     Compute the singular values of (I - Q Q') A
      Call dgesvd('N','N',m,n,imqqa,lda,s,u,ldu,vt,ldvt,work,lwork,info)

!     Calculate probabilistic error estimate
      Allocate (omega(n),y(m),res(m),res_norm(r))

      Do j = 1, r
!       Generate Gaussian random vector, omega
        Call g05skf(n,zero,one,state,omega,ifail)
!       calculate y = A omega
        Call dgemv('N',m,n,one,a,m,omega,1,zero,y,1)
!       Calculate (I - Q Q') y
        Call dgemv('N',m,m,one,imqq,m,y,1,zero,res,1)
        res_norm(j) = dnrm2(m,res,1)
      End Do
      pi = x01aaf(pi)
      err_bound = sqrt(200.0_nag_wp/pi)*maxval(res_norm)

!     Largest singular value = 2-norm of (I - Q Q') A
      Write (nout,*)
      Write (nout,*) 'True error'
      Write (nout,99998) s(1)
      Write (nout,*)
!     Error bound
      Write (nout,*) 'Probabilistic error bound'
      Write (nout,99998) err_bound
      Write (nout,*) 'Maximum probability that true error exceeds error bound'
      Write (nout,99998)(10.0_nag_wp**(-r))

      Deallocate (state)
      Deallocate (work)
      Deallocate (a)
      Deallocate (pa)
      Deallocate (imqq)
      Deallocate (imqqa)
      Deallocate (res)
      Deallocate (res_norm)
      Deallocate (tau)
      Deallocate (q)
      Deallocate (s)
      Deallocate (omega)
      Deallocate (y)

99998 Format (1X,(3X,E11.4))
    End Program f10dafe