Program c09ecfe
! C09ECF Example Program Text
! Mark 28.4 Release. NAG Copyright 2022.
! .. 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, i1, ifail, ilevel, itype_coeffs, &
j1, lda, ldb, ldcoefs, lenc, m, n, &
nf, nwcn, nwct, nwl, nwlinv, nwlmax
Character (10) :: mode, wavnam, wtrans
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: a(:,:), b(:,:), c(:), coefs(:,:)
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,nwlmax,nf,nwct,nwcn,icomm,ifail)
lenc = nwct
Allocate (c(lenc),dwtlvm(nwlmax),dwtlvn(nwlmax))
nwl = nwlmax
! 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)
! Allocate an array in which to extract coefficients
ldcoefs = dwtlvm(nwl)
Allocate (coefs(ldcoefs,dwtlvn(nwl)))
! Extract each set of coefficients, working from the deepest level
Write (nout,99993)
Do ilevel = nwl, 1, -1
Write (nout,99992) ilevel, dwtlvm(nwl-ilevel+1), dwtlvn(nwl-ilevel+1)
Do itype_coeffs = 0, 3
Select Case (itype_coeffs)
Case (0)
If (ilevel==nwl) Then
Write (nout,99991) 'Approximation coefficients '
End If
Case (1)
Write (nout,99991) 'Vertical coefficients '
Case (2)
Write (nout,99991) 'Horizontal coefficients '
Case (3)
Write (nout,99991) 'Diagonal coefficients '
End Select
If (itype_coeffs>0 .Or. ilevel==nwl) Then
! Call the 2D extraction routine c09eaf
Call c09eyf(ilevel,itype_coeffs,lenc,c,coefs,ldcoefs,icomm,ifail)
Do i1 = 1, dwtlvm(nwl-ilevel+1)
Write (nout,99989)(coefs(i1,j1),j1=1,dwtlvn(nwl-ilevel+1))
End Do
End If
End Do
End Do
nwlinv = nwl
! Reconstruct original data
ifail = 0
Call c09edf(nwlinv,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