Program f01sbfe
! F01SBF Example Program Text
! Mark 30.2 Release. NAG Copyright 2024.
! .. Use Statements ..
Use nag_library, Only: dgemm, f01sbf, f06raf, f11mkf, nag_wp, x04caf
! .. 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 :: nin = 5, nout = 6
! .. Local Scalars ..
Real (Kind=nag_wp) :: errtol, norm, norma
Integer :: element, i, icount, ifail, irevcm, &
j, k, lda_dense, ldh, ldht, ldw, m, &
maxit, n, nnz, q, seed
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: a(:), a_dense(:,:), comm(:), h(:,:), &
ht(:,:), w(:,:)
Real (Kind=nag_wp) :: work(1)
Integer, Allocatable :: icolzp(:), irowix(:)
Integer :: icomm(9)
! .. Intrinsic Procedures ..
Intrinsic :: max
! .. Executable Statements ..
Write (nout,*) 'F01SBF Example Program Results'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
Read (nin,*) m, n, k
! Pretend A is a square q x q matrix to allow use of f11mkf for matrix multiplications
q = max(m,n)
ldw = q
ldh = k
ldht = q
lda_dense = m
Allocate (icolzp(q+1))
Allocate (w(ldw,n))
Allocate (h(ldh,n))
Allocate (ht(ldht,k))
Allocate (comm((2*m+n)*k+3))
Allocate (a_dense(lda_dense,n))
! Initialize W and Ht to zero
w(1:ldw,1:n) = zero
ht(1:ldht,1:k) = zero
! Read A from data file
Read (nin,*) icolzp(1:n+1)
nnz = icolzp(n+1) - 1
Allocate (a(nnz),irowix(nnz))
Read (nin,*)(a(i),irowix(i),i=1,nnz)
! If m > n (q==m), f11mkf will need a total of q columns
! otherwise no adjustment is necessary
icolzp(n+2:q+1) = icolzp(n+1)
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = -1
icount = 0
irevcm = 0
! Choose values for the arguments errtol and seed
errtol = 1.0E-3_nag_wp
seed = 234
! We will terminate after 200 iterations
maxit = 200
! Find a non-negative matrix factorixation A ~= WH
revcomm: Do
Call f01sbf(irevcm,m,n,k,w,ldw,h,ldh,ht,ldht,seed,errtol,comm,icomm, &
ifail)
Select Case (irevcm)
Case (1)
If (icount==maxit) Then
Write (nout,*) 'Exiting after the maximum number of iterations'
Exit revcomm
End If
icount = icount + 1
! Insert print statements here to monitor progress if desired.
Cycle revcomm
Case (2)
! Compute a^t * w and store in ht
Call f11mkf('T',q,k,one,icolzp,irowix,a,w,ldw,zero,ht,ldht,ifail)
Cycle revcomm
Case (3)
! Compute a * ht and store in w
Call f11mkf('N',q,k,one,icolzp,irowix,a,ht,ldht,zero,w,ldw,ifail)
Cycle revcomm
Case Default
Exit revcomm
End Select
End Do revcomm
If (ifail==0 .Or. icount==maxit) Then
! Print solution
Write (nout,*)
Call x04caf('General',' ',m,k,w,ldw,'W',ifail)
Write (nout,*)
Call x04caf('General',' ',k,n,h,ldh,'H',ifail)
! In order to compute the residual we'll convert a to dense format
a_dense(1:lda_dense,1:n) = zero
element = 1
Do i = 1, n
Do j = 1, icolzp(i+1) - icolzp(i)
a_dense(irowix(element),i) = a(element)
element = element + 1
End Do
End Do
norma = f06raf('F',m,n,a_dense,lda_dense,work)
Call dgemm('n','n',m,n,k,-one,w,ldw,h,ldh,one,a_dense,lda_dense)
Write (nout,*)
norm = f06raf('F',m,n,a_dense,lda_dense,work)
Write (nout,*)
Write (nout,99999) 'The relative residual norm, ||A-WH||/||A||, is: ', &
norm/norma
End If
99999 Format (1X,A,Es9.2)
End Program f01sbfe