Program f10cafe

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

!     .. Use Statements ..
      Use nag_library, Only: dgesvd, f10caf, g05kff, nag_wp, x02ajf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: lseed = 1, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: rtol_abs, rtol_rel
      Integer                          :: genid, i, ifail, info, k, lda, ldu,  &
                                          ldvt, lstate, lwork, m, n, r, subid
      Character (1)                    :: jobu, jobvt
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), a_copy(:,:), s(:), u(:,:),   &
                                          vt(:,:), work(:)
      Real (Kind=nag_wp)               :: dummy(1)
      Integer                          :: seed(lseed)
      Integer, Allocatable             :: state(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: min, nint
!     .. Executable Statements ..
      Write (nout,*) 'F10CAF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) m, n
      lda = m
      ldu = m
      ldvt = n
      Allocate (a(lda,n),a_copy(m,n),s(n),u(ldu,m))

!     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 arguments to randmoised SVD routine (f10caf) - jobu, jobvt, k
      Read (nin,*) jobu, jobvt, k

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

!     Evaluate SVD using deterministic algorithm

      Allocate (vt(n,n))
!     Use routine workspace query to get optimal workspace.
      lwork = -1
!     The NAG name equivalent of dgesvd is f08kbf
      Call dgesvd('A','S',m,n,a,lda,s,u,ldu,vt,ldvt,dummy,lwork,info)

!     Make sure that there is enough workspace
      lwork = nint(dummy(1))
      Allocate (work(lwork))

!     Compute the singular values and left and right singular vectors
!     of A.

!     The NAG name equivalent of dgesvd is f08kbf
      Call dgesvd('A','S',m,n,a,lda,s,u,ldu,vt,ldvt,work,lwork,info)

      If (info/=0) Then
        Write (nout,99998) 'Failure in DGESVD. INFO =', info
99998   Format (1X,A,I4)
        Go To 100
      End If

!     Print the significant singular values of A
      Write (nout,*)
      Write (nout,*) 'Singular values from deterministic SVD (dgesvd):'
      Write (nout,99997) s(1:min(m,n))

99997 Format (6(1X,F9.4))

!     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 randomised SVD
      rtol_abs = x02ajf()**0.875_nag_wp
      rtol_rel = rtol_abs

      ifail = 0
      Call f10caf(jobu,jobvt,m,n,a_copy,lda,k,rtol_abs,rtol_rel,state,s,u,ldu, &
        vt,ldvt,r,ifail)

      Write (nout,*)
      Write (nout,99996) 'Rank estimate from randomized SVD (f10caf) = ', r
      Write (nout,*)
      Write (nout,*) 'Singular values from randomized SVD'
      Write (nout,99997) s(1:r)
99996 Format (1X,A,I0)
100   Continue
    End Program f10cafe