Program f10dafe
! F10DAF Example Program
! Mark 28.5 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