PROGRAM c09ecfe

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : c09abf, c09ecf, c09edf, nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: nin = 5, nout = 6
!      .. Local Scalars ..
       INTEGER                         :: i, i1, i2, ifail, ilevel, iskip,     &
                                          itype_coeffs, j1, jstart, lda, ldb,  &
                                          lenc, m, n, nf, nwcn, nwct, nwl
       CHARACTER (10)                  :: mode, wavnam, wtrans
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: a(:,:), b(:,:), c(:)
       INTEGER, ALLOCATABLE            :: dwtlvm(:), dwtlvn(:)
       INTEGER                         :: icomm(180)
!      .. Executable Statements ..
       WRITE (nout,*) 'C09ECF Example Program Results'
       WRITE (nout,*)
!      Skip heading in data file
       READ (nin,*)
!      Read problem parameters
       READ (nin,*) m, n
       READ (nin,*) wavnam, mode
       lda = m
       ldb = m
       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 = 0
       CALL c09abf(wavnam,wtrans,mode,m,n,nwl,nf,nwct,nwcn,icomm,ifail)

       lenc = nwct
       ALLOCATE (c(lenc),dwtlvm(nwl),dwtlvn(nwl))

!      Perform Discrete Wavelet transform
       ifail = 0
       CALL c09ecf(m,n,a,lda,lenc,c,nwl,dwtlvm,dwtlvn,icomm,ifail)

       WRITE (nout,99997) nwl
       WRITE (nout,99996)
       WRITE (nout,99995) dwtlvm(1:nwl)
       WRITE (nout,99994)
       WRITE (nout,99995) dwtlvn(1:nwl)

       WRITE (nout,99993)
       jstart = 1
       DO ilevel = nwl, 1, -1
          WRITE (nout,99992) ilevel, dwtlvm(nwl-ilevel+1), &
             dwtlvn(nwl-ilevel+1)

          iskip = dwtlvm(nwl-ilevel+1)
          i2 = iskip*dwtlvn(nwl-ilevel+1) - 1

          DO itype_coeffs = 1, 4
             SELECT CASE (itype_coeffs)
             CASE (1)
                IF (ilevel==nwl) THEN
                   WRITE (nout,99991) 'Approximation coefficients '
                END IF
             CASE (2)
                WRITE (nout,99991) 'Vertical coefficients      '
             CASE (3)
                WRITE (nout,99991) 'Horizontal coefficients    '
             CASE (4)
                WRITE (nout,99991) 'Diagonal coefficients      '
             END SELECT
             IF (itype_coeffs>1 .OR. ilevel==nwl) THEN
                DO i1 = jstart, jstart + iskip - 1
                   WRITE (nout,99989) (c(j1),j1=i1,i1+i2,iskip)
                END DO
                jstart = jstart + i2 + 1
             END IF
          END DO
       END DO

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

       WRITE (nout,99990)
       DO i = 1, m
          WRITE (nout,99998) b(i,1:n)
       END DO

99999  FORMAT (1X,' MLDWT :: Wavelet  : ',A/1X,'          End mode : ',A/1X, &
          '          M        : ',I4/1X,'          N        : ',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,' Wavelet coefficients C : ')
99992  FORMAT (1X,55('-')/1X,' Level : ',I10,'; output is ',I10,' by ',I10/1X, &
          55('-'))
99991  FORMAT (1X,A28,': ')
99990  FORMAT (/1X,' Reconstruction           B : ')
99989  FORMAT (4X,5(F8.4,1X):)
    END PROGRAM c09ecfe