Program c09abfe

!     C09ABF Example Program Text
!     Mark 26.1 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: c09abf, c09ecf, c09edf, c09eyf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Integer                          :: i, ifail, lda, ldb, lenc, m, n, nf,  &
                                          nwcm, nwcn, nwct, nwl, nwlmax,       &
                                          want_coeffs, want_level
      Character (10)                   :: mode, wavnam, wtrans
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), b(:,:), c(:), d(:,:)
      Integer, Allocatable             :: dwtlevm(:), dwtlevn(:)
      Integer                          :: icomm(180)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: sum
!     .. Executable Statements ..
      Continue
      Write (nout,*) 'C09ABF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
!     Read problem parameters
      Read (nin,*) m, n
      lda = m
      ldb = m
      Read (nin,*) wavnam, mode
      Allocate (a(lda,n),b(ldb,n))

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

!     Read data array and write it out
      Do i = 1, m
        Read (nin,*) a(i,1:n)
      End Do

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

!     Query wavelet filter dimensions
!     For Multi-Resolution Analysis, decomposition, wtrans = 'M'
      wtrans = 'Multilevel'

!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call c09abf(wavnam,wtrans,mode,m,n,nwlmax,nf,nwct,nwcn,icomm,ifail)

      lenc = nwct
      Allocate (c(lenc),dwtlevm(nwlmax),dwtlevn(nwlmax))

!     Calculate one less than the max possible number of levels
      nwl = nwlmax - 1

!     Perform Discrete Wavelet transform
      ifail = 0
      Call c09ecf(m,n,a,lda,lenc,c,nwl,dwtlevm,dwtlevn,icomm,ifail)

!     c09abf returns nwct based on max levels, so recalculate.
      nwct = sum(3*dwtlevm(1:nwl)*dwtlevn(1:nwl))
      nwct = nwct + dwtlevm(1)*dwtlevn(1)

!     Print the details of the transform.
      Write (nout,*)
      Write (nout,99997) nf
      Write (nout,99996) nwl
      Write (nout,99995)
      Write (nout,99994) dwtlevm(1:nwl)
      Write (nout,99993)
      Write (nout,99994) dwtlevn(1:nwl)
      Write (nout,99992) nwct
      Write (nout,*)
      Write (nout,99991)
      Write (nout,99998) c(1:nwct)

!     Now select a nominated matrix of coefficients at a nominated level.
!     Remember that level 0 is input data, 1 first coeffs and so on up to nwl,
!     which is the deepest level and contains approx. coefficients.
      want_level = nwl - 1
!     Print only vertical detail coeffs at selected level.
      want_coeffs = 1
      nwcm = dwtlevm(nwl-want_level+1)
      nwcn = dwtlevn(nwl-want_level+1)
      Allocate (d(nwcm,nwcn))

!     Extract the selected set of coefficients.
      Call c09eyf(want_level,want_coeffs,lenc,c,d,nwcm,icomm,ifail)

!     Print out the selected coefficients
      Write (nout,*)
      Write (nout,99989) want_coeffs, want_level
      Do i = 1, nwcm
        Write (nout,99998) d(i,1:nwcn)
      End Do

!     Reconstruct original data
      ifail = 0
      Call c09edf(nwl,lenc,c,m,n,b,ldb,icomm,ifail)

      Write (nout,*)
      Write (nout,99990)
      Do i = 1, m
        Write (nout,99998) b(i,1:n)
      End Do

99999 Format (1X,' Parameters read from file :: ',/,'     Wavelet : ',A10,     &
        ' End mode : ',A10,' M = ',I10,' N = ',I10)
99998 Format (8(F8.4,1X),:)
99997 Format (1X,' Length of wavelet filter :             ',I10)
99996 Format (1X,' Number of Levels :                     ',I10)
99995 Format (1X,                                                              &
        ' Number of coefficients in first dimension for each level : ')
99994 Format (16X,8(I8,1X),:)
99993 Format (1X,                                                              &
        ' Number of coefficients in second dimension for each level : ')
99992 Format (1X,' Total number of wavelet coefficients : ',I10)
99991 Format (1X,' All Wavelet coefficients C : ')
99990 Format (1X,' Reconstruction       B : ')
99989 Format (1X,' Type ',I1,' coefficients at selected wavelet level, ',I4,   &
        ': ')
    End Program c09abfe