    Program c09acfe

!     Mark 24 Release. NAG Copyright 2012.

!     .. Use Statements ..
      Use nag_library, Only: c09acf, c09fcf, c09fdf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Integer                          :: fr, i, i1, ifail, ilevel,            &
                                          itype_coeffs, j, k, lda, ldb, ldd,   &
                                          lenc, lmax, locc, m, n, nf, nwcfr,   &
                                          nwcm, nwcn, nwct, nwl, sda, sdb,     &
                                          sdd, want_coeffs, want_level
      Character (10)                   :: mode, wavnam, wtrans
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:,:), b(:,:,:), c(:), d(:,:,:)
      Integer, Allocatable             :: dwtlvfr(:), dwtlvm(:), dwtlvn(:)
      Integer                          :: icomm(260)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: sum
!     .. Executable Statements ..
      Continue
      Write (nout,*) 'C09ACF 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))

      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) Read (nin,*)
      End Do

      Write (nout,*) ' Input Data     A :'
      Do j = 1, fr
        Do i = 1, m
          Write (nout,99998) a(i,1:n,j)
        End Do
        Write (nout,*)
      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,lmax,nf,nwct,nwcn,nwcfr,icomm, &
        ifail)

!     Transform one less than the max possible number of levels.
      nwl = lmax - 1

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

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

!     c09acf returns nwct based on max levels, so recalculate.
      nwct = sum(7*dwtlvm(1:nwl)*dwtlvn(1:nwl)*dwtlvfr(1:nwl))
      nwct = nwct + dwtlvm(1)*dwtlvn(1)*dwtlvfr(1)

      Write (nout,99997) nwl
      Write (nout,99988) nf
      Write (nout,99987) nwct
      Write (nout,99996)
      Write (nout,99995) dwtlvm(1:nwl)
      Write (nout,99994)
      Write (nout,99995) dwtlvn(1:nwl)
      Write (nout,99993)
      Write (nout,99995) dwtlvfr(1:nwl)

!     Select the deepest level.
      want_level = nwl
!     Select the approximation coefficients.
      want_coeffs = 0

!     Identify each set of coefficients in C
      Do ilevel = nwl, 1, -1

        If (ilevel/=want_level) Cycle

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

        ldd = nwcm
        sdd = nwcn
        Allocate (d(ldd,sdd,nwcfr))

        Write (nout,99986) ilevel, nwcm, nwcn, nwcfr

        Do itype_coeffs = 0, 7

          If (itype_coeffs/=want_coeffs) Cycle

!         Unless we're looking at the deepest level of nesting, which contains
!         approximation coefficients, advance the pointer on past the preceding
!         levels
          If (ilevel==nwl) Then
            locc = 0
          Else
            locc = 8*dwtlvm(1)*dwtlvn(1)*dwtlvfr(1)
            Do i = ilevel + 1, nwl - 1
              locc = locc + 7*dwtlvm(nwl-i+1)*dwtlvn(nwl-i+1)*dwtlvfr(nwl-i+1)
            End Do
          End If

!         Now decide which coefficient type we are considering
          Select Case (itype_coeffs)
          Case (0)
            If (ilevel==nwl) Then
              Write (nout,99985) 'Approximation coefficients (LLL) '
              locc = locc + 1
            End If
          Case (1)
            Write (nout,99985) 'Detail coefficients (LLH)  '
            If (ilevel==nwl) Then
!             Advance pointer past approximation coefficients
              locc = locc + nwcm*nwcn*nwcfr + 1
            Else
              locc = locc + 1
            End If
          Case (2)
            Write (nout,99985) 'Detail coefficients (LHL)  '
            If (ilevel==nwl) Then
!             Advance pointer past approximation coefficients and 1 set of
!             detail coefficients
              locc = locc + 2*nwcm*nwcn*nwcfr + 1
            Else
!             Advance pointer past 1 set of detail coefficients
              locc = locc + nwcm*nwcn*nwcfr + 1
            End If
          Case (3)
            Write (nout,99985) 'Detail coefficients (LHH)  '
            If (ilevel==nwl) Then
!             Advance pointer past approximation coefficients and 2 sets of
!             detail coefficients
              locc = locc + 3*nwcm*nwcn*nwcfr + 1
            Else
!             Advance pointer past 2 sets of detail coefficients
              locc = locc + 2*nwcm*nwcn*nwcfr + 1
            End If
          Case (4)
            Write (nout,99985) 'Detail coefficients (HLL)  '
            If (ilevel==nwl) Then
!             Advance pointer past approximation coefficients and 3 sets of
!             detail coefficients
              locc = locc + 4*nwcm*nwcn*nwcfr + 1
            Else
!             Advance pointer past 3 sets of detail coefficients
              locc = locc + 3*nwcm*nwcn*nwcfr + 1
            End If
          Case (5)
            Write (nout,99985) 'Detail coefficients (HLH)  '
            If (ilevel==nwl) Then
!             Advance pointer past approximation coefficients and 4 sets of
!             detail coefficients
              locc = locc + 5*nwcm*nwcn*nwcfr + 1
            Else
!             Advance pointer past 4 sets of detail coefficients
              locc = locc + 4*nwcm*nwcn*nwcfr + 1
            End If
          Case (6)
            Write (nout,99985) 'Detail coefficients (HHL)  '
            If (ilevel==nwl) Then
!             Advance pointer past approximation coefficients and 5 sets of
!             detail coefficients
              locc = locc + 6*nwcm*nwcn*nwcfr + 1
            Else
!             Advance pointer past 4 sets of detail coefficients
              locc = locc + 5*nwcm*nwcn*nwcfr + 1
            End If
          Case (7)
            Write (nout,99985) 'Detail coefficients (HHH)  '
            If (ilevel==nwl) Then
!             Advance pointer past approximation coefficients and 6 sets of
!             detail coefficients
              locc = locc + 7*nwcm*nwcn*nwcfr + 1
            Else
!             Advance pointer past 5 sets of detail coefficients
              locc = locc + 6*nwcm*nwcn*nwcfr + 1
            End If
          End Select

          If (itype_coeffs>0 .Or. ilevel==nwl) Then

            If (itype_coeffs==0) Then

!             For a multi level transform approx coeffs stored as nwcm x nwcn x nwcfr...
              i1 = locc
              Do k = 1, nwcfr
                Do j = 1, nwcn
                  Do i = 1, nwcm
                    d(i,j,k) = c(i1)
                    i1 = i1 + 1
                  End Do
                End Do
              End Do

            Else

!             but detail coefficients are stored as ncwfr x nwcm x nwcn
              Do k = 1, nwcfr
                Do j = 1, nwcn
                  Do i = 1, nwcm
                    i1 = locc - 1 + (j-1)*nwcfr*nwcm + (i-1)*nwcfr + k
                    d(i,j,k) = c(i1)
                  End Do
                End Do
              End Do

            End If

!           Print out the selected set of coefficients
            Write (nout,99989) ilevel, itype_coeffs
            Do k = 1, nwcfr
              Write (nout,99990) k
              Do i = 1, nwcm
                Write (nout,99991) d(i,1:nwcn,k)
              End Do
            End Do

          End If
        End Do
        Deallocate (d)
      End Do

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

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

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,' Number of Levels :         ',I10)
99996 Format (1X,' Number of coefficients in 1st dimension for each level :')
99995 Format (8(I8,1X):)
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 (/1X,' Reconstruction           B : ')
99991 Format (1X,8(F8.4,1X):)
99990 Format (1X,' Frame ',I2,' : ')
99989 Format (1X,' Level ',I2,', Coefficients ',I2,' : ')
99988 Format (1X,' Length of wavelet filter : ',I10)
99987 Format (1X,' Total number of wavelet coefficients : ',I10)
99986 Format (/1X,70('-')/1X,'Level : ',I10,'; output is ',I10,' by ',I10, &
        ' by ',I10/1X,70('-'))
99985 Format (/1X,A)
    End Program c09acfe
