NAG Library Manual, Mark 28.7
```    Program f01sbfe

!     F01SBF Example Program Text

!     Mark 28.7 Release. NAG Copyright 2022.

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

!     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
nnz = icolzp(n+1) - 1
Allocate (a(nnz),irowix(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
```