Program c09fcfe

!     Mark 24 Release. NAG Copyright 2012.

!     .. Use Statements ..
      Use nag_library, Only: c09acf, c09fcf, c09fdf, 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, i1, ifail, ilevel,            &
                                          itype_coeffs, j, k, lda, ldb, ldd,   &
                                          lenc, 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(:,:,:),  &
                                          e(:,:,:)
      Integer, Allocatable             :: dwtlvfr(:), dwtlvm(:), dwtlvn(:)
      Integer                          :: icomm(260)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max, real, sqrt
!     .. Executable Statements ..
      Continue
      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) Read (nin,*)
      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,nwl,nf,nwct,nwcn,nwcfr,icomm, &
        ifail)

      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)

      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

!     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,99987) 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,99986) 'Approximation coefficients (LLL) '
              locc = locc + 1
            End If
          Case (1)
            Write (nout,99986) '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,99986) '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,99986) '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,99986) '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,99986) '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,99986) '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,99986) '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,99988) k
              Do i = 1, nwcm
                Write (nout,99998) 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)

!     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