Program f01sbfe

!     F01SBF Example Program Text

!     Mark 28.4 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
      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