Program c09fcfe

!     C09FCF Example Program Text
!     Mark 28.5 Release. NAG Copyright 2022.

!     .. Use Statements ..
      Use nag_library, Only: c09acf, c09fcf, c09fdf, c09fyf, nag_wp, x02ajf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: eps, esq, frob
      Integer                          :: fr, i, ifail, j, k, lda, ldb, ldd,   &
                                          lenc, m, n, nf, nwcfr, nwcm, nwcn,   &
                                          nwct, nwl, nwlinv, nwlmax, sda, sdb, &
                                          sdd, want_coeffs, want_level
      Character (10)                   :: mode, wavnam, wtrans
      Character (33)                   :: title
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:,:), b(:,:,:), c(:), d(:,:,:),  &
                                          e(:,:,:)
      Integer, Allocatable             :: dwtlvfr(:), dwtlvm(:), dwtlvn(:)
      Integer                          :: icomm(260)
      Character (3)                    :: cpass(0:7)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max, real, sqrt
!     .. Executable Statements ..
      Write (nout,*) 'C09FCF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
!     Read problem parameters
      Read (nin,*) m, n, fr
      Read (nin,*) wavnam, mode
      lda = m
      sda = n
      ldb = m
      sdb = n
      Allocate (a(lda,sda,fr),b(ldb,sdb,fr),e(m,n,fr))

      Write (nout,99999) wavnam, mode, m, n, fr

!     Read data array and write it out

      Do j = 1, fr
        Do i = 1, m
          Read (nin,*) a(i,1:n,j)
        End Do
        If (j<fr) Then
          Read (nin,*)
        End If
      End Do

      Write (nout,*) ' Input Data     A :'
      Do j = 1, fr
        Write (nout,99997) j
        Do i = 1, m
          Write (nout,99998) a(i,1:n,j)
        End Do
      End Do

!     Query wavelet filter dimensions
!     For Multi-Resolution Analysis, decomposition, wtrans = 'M'
      wtrans = 'Multilevel'
      ifail = 0
      Call c09acf(wavnam,wtrans,mode,m,n,fr,nwlmax,nf,nwct,nwcn,nwcfr,icomm,   &
        ifail)

      lenc = nwct
      Allocate (c(lenc),dwtlvm(nwlmax),dwtlvn(nwlmax),dwtlvfr(nwlmax))

      nwl = nwlmax

!     Perform Discrete Wavelet transform
      ifail = 0
      Call c09fcf(m,n,fr,a,lda,sda,lenc,c,nwl,dwtlvm,dwtlvn,dwtlvfr,icomm,     &
        ifail)

      Write (nout,99996) nwl
      Write (nout,99995)
      Write (nout,99992) dwtlvm(1:nwl)
      Write (nout,99994)
      Write (nout,99992) dwtlvn(1:nwl)
      Write (nout,99993)
      Write (nout,99992) dwtlvfr(1:nwl)

!     Print the first level HLL coefficients
      want_level = 1
      want_coeffs = 4

      nwcm = dwtlvm(nwl-want_level+1)
      nwcn = dwtlvn(nwl-want_level+1)
      nwcfr = dwtlvfr(nwl-want_level+1)

!     Allocate space to store the selected coefficients
      ldd = nwcm
      sdd = nwcn
      Allocate (d(ldd,sdd,nwcfr))

      Write (nout,99987) want_level, nwcm, nwcn, nwcfr

      cpass(0:7) = (/'LLL','LLH','LHL','LHH','HLL','HLH','HHL','HHH'/)
      If (want_coeffs==0) Then
        title = 'Approximation coefficients (LLL)'
      Else
        title = 'Detail coefficients (' // cpass(want_coeffs) // ')'
      End If

!     Extract coefficients
      Call c09fyf(want_level,want_coeffs,lenc,c,d,ldd,sdd,icomm,ifail)

!     Print out the selected set of coefficients
      Write (nout,99986) title
      Write (nout,99989) want_level, want_coeffs
      Do k = 1, nwcfr
        Write (nout,99988) k
        Do i = 1, nwcm
          Write (nout,99998) d(i,1:nwcn,k)
        End Do
      End Do

      nwlinv = nwl

!     Reconstruct original data
      ifail = 0
      Call c09fdf(nwlinv,lenc,c,m,n,fr,b,ldb,sdb,icomm,ifail)

!     Check reconstruction matches original
      eps = 10.0_nag_wp*real(m,kind=nag_wp)*real(n,kind=nag_wp)*               &
        real(fr,kind=nag_wp)*x02ajf()

      e(1:m,1:n,1:fr) = b(1:m,1:n,1:fr) - a(1:m,1:n,1:fr)
      frob = 0.0_nag_wp
      Do k = 1, fr
        esq = 0.0_nag_wp
        Do j = 1, n
          Do i = 1, m
            esq = esq + e(i,j,k)**2
          End Do
        End Do
        frob = max(frob,sqrt(esq))
      End Do

      If (frob>eps) Then
        Write (nout,99991)
      Else
        Write (nout,99990)
      End If

99999 Format (1X,' MLDWT :: Wavelet  : ',A,/,1X,'          End mode : ',A,/,   &
        1X,'          M        : ',I4,/,1X,'          N        : ',I4,/,1X,    &
        '          FR       : ',I4,/)
99998 Format (8(F8.4,1X),:)
99997 Format (1X,' Frame ',I2,' : ')
99996 Format (/,1X,' Number of Levels : ',I10)
99995 Format (1X,' Number of coefficients in 1st dimension for each level :')
99994 Format (1X,' Number of coefficients in 2nd dimension for each level :')
99993 Format (1X,' Number of coefficients in 3rd dimension for each level :')
99992 Format (8(I8,1X),:)
99991 Format (/,1X,' Fail: Frobenius norm of B-A, where A is the original ',/, &
        1X,' data and B is the reconstrucion, is too large.')
99990 Format (/,1X,' Success: the reconstruction matches the original.')
99989 Format (1X,' Level ',I2,', Coefficients ',I2,' : ')
99988 Format (1X,' Frame ',I2,' : ')
99987 Format (/,1X,70('-'),/,1X,'Level : ',I10,'; output is ',I10,' by ',I10,  &
        ' by ',I10,/,1X,70('-'))
99986 Format (/,1X,A)
    End Program c09fcfe