Program c09fafe

!     Mark 24 Release. NAG Copyright 2012.

!     .. Use Statements ..
      Use nag_library, Only: c09acf, c09faf, c09fbf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Integer                          :: cindex, fr, i, i1, ifail, j, k, lda, &
                                          ldb, lenc, locc, m, n, nf, nwcfr,    &
                                          nwcm, nwcn, nwct, nwl, sda, sdb
      Character (12)                   :: mode, wavnam, wtrans
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:,:), b(:,:,:), c(:), d(:,:,:)
      Integer                          :: icomm(260)
!     .. Executable Statements ..
      Continue
      Write (nout,*) 'C09FAF Example Program Results'

!     Skip heading in data file
      Read (nin,*)
!     Read problem parameters.
      Read (nin,*) m, n, fr
      Read (nin,*) wavnam, mode
      Write (nout,99999) wavnam, mode

      lda = m
      sda = n
      Allocate (a(lda,sda,fr))
      ldb = m
      sdb = n
      Allocate (b(ldb,sdb,fr))

!     Read data array
      Do j = 1, fr
        Do i = 1, m
          Read (nin,*) a(i,1:n,j)
        End Do
        Read (nin,*)
      End Do

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

!     Query wavelet filter dimensions
      wtrans = 'Single Level'

!     ifail: behaviour on error exit   
!     =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call c09acf(wavnam,wtrans,mode,m,n,fr,nwl,nf,nwct,nwcn,nwcfr,icomm, &
        ifail)
      nwcm = nwct/(8*nwcn*nwcfr)
      lenc = nwct
      Allocate (c(lenc))

!     3D DWT decomposition
      ifail = 0
      Call c09faf(m,n,fr,a,lda,sda,lenc,c,icomm,ifail)

      Allocate (d(nwcm,nwcn,nwcfr))

      Do cindex = 0, 7

!         Decide which coefficient type we are considering
!         and advance the pointer locc to the first element 
!         of that 3D array in C.
        Select Case (cindex)
        Case (0)
          Write (nout,99990) 'Approximation coefficients (LLL) '
          locc = 1
        Case (1)
          Write (nout,99990) 'Detail coefficients (LLH)  '
!         Advance pointer past approximation coefficients
          locc = nwcm*nwcn*nwcfr + 1
        Case (2)
          Write (nout,99990) 'Detail coefficients (LHL)  '
!             Advance pointer past approximation coefficients and 1 set of
!             detail coefficients
          locc = 2*nwcm*nwcn*nwcfr + 1
        Case (3)
          Write (nout,99990) 'Detail coefficients (LHH)  '
!             Advance pointer past approximation coefficients and 2 sets of
!             detail coefficients
          locc = 3*nwcm*nwcn*nwcfr + 1
        Case (4)
          Write (nout,99990) 'Detail coefficients (HLL)  '
!             Advance pointer past approximation coefficients and 3 sets of
!             detail coefficients
          locc = 4*nwcm*nwcn*nwcfr + 1
        Case (5)
          Write (nout,99990) 'Detail coefficients (HLH)  '
!             Advance pointer past approximation coefficients and 4 sets of
!             detail coefficients
          locc = 5*nwcm*nwcn*nwcfr + 1
        Case (6)
          Write (nout,99990) 'Detail coefficients (HHL)  '
!             Advance pointer past approximation coefficients and 5 sets of
!             detail coefficients
          locc = 6*nwcm*nwcn*nwcfr + 1
        Case (7)
          Write (nout,99990) 'Detail coefficients (HHH)  '
!             Advance pointer past approximation coefficients and 6 sets of
!             detail coefficients
          locc = 7*nwcm*nwcn*nwcfr + 1
        End Select

        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

        Do j = 1, nwcfr
          If (j==1) Then
            Write (nout,99995,Advance='no') j
          Else
            Write (nout,99994,Advance='no') j
          End If
        End Do
        Write (nout,*)
        Do i = 1, nwcm
          Do j = 1, nwcfr
            If (i==1 .And. j==1) Then
              Write (nout,99993,Advance='no') cindex, d(i,1:nwcn,j)
            Else If (j==1) Then
              Write (nout,99992,Advance='no') d(i,1:nwcn,j)
            Else
              Write (nout,99991,Advance='no') d(i,1:nwcn,j)
            End If
          End Do
          Write (nout,*)
        End Do
      End Do

!     3D DWT reconstruction
      ifail = 0
      Call c09fbf(m,n,fr,lenc,c,b,ldb,sdb,icomm,ifail)

      Write (nout,99998) 'Output Data                    B'
      Do j = 1, fr
        Write (nout,99996) j
        Do i = 1, m
          Write (nout,99997) b(i,1:n,j)
        End Do
      End Do

99999 Format (/1X,'DWT ::'/1X,'       Wavelet : ',A/1X,'       End mode: ',A)
99998 Format (/1X,A,' : ')
99997 Format (1X,8(F8.4,1X):)
99996 Format (1X,'Frame ',I2,' : ')
99995 Format (1X,'Coefficients        Frame ',I2:)
99994 Format (14X,'Frame ',I2:)
99993 Format (4X,I4,7X,8(F8.4,1X):)
99992 Format (15X,8(F8.4,1X):)
99991 Format (1X,3X,8(F8.4,1X):)
99990 Format (/1X,A)
    End Program c09fafe