Program f01safe

!     F01SAF Example Program Text

!     Mark 28.4 Release. NAG Copyright 2022.

!     .. Use Statements ..
      Use nag_library, Only: dgemm, f01saf, f06raf, nag_wp, x04caf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: errtol, norm, norma
      Integer                          :: i, ifail, k, lda, ldh, ldw, m,       &
                                          maxit, n, seed
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), h(:,:), w(:,:)
      Real (Kind=nag_wp)               :: work(1)
!     .. Executable Statements ..
      Write (nout,*) 'F01SAF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) m, n, k
      lda = m
      ldw = m
      ldh = k
      Allocate (a(lda,n))
      Allocate (w(ldw,n))
      Allocate (h(ldh,n))
!     Read A from data file
      Read (nin,*)(a(i,1:n),i=1,m)

!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = -1

!     Choose values for errtol and seed
      errtol = 1.0E-6_nag_wp
      seed = 23

!     Use the default value for maxit
      maxit = -1

!     Find the non-negative matrix factorixation A ~= WH
      Call f01saf(m,n,k,a,lda,w,ldw,h,ldh,seed,errtol,maxit,ifail)

      If (ifail==0) 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)

        norma = f06raf('F',m,n,a,lda,work)

        Call dgemm('n','n',m,n,k,-1.0_nag_wp,w,ldw,h,ldh,1.0_nag_wp,a,lda)
        Write (nout,*)
        norm = f06raf('F',m,n,a,lda,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 f01safe