Program c09abfe
! Mark 24 Release. NAG Copyright 2012.
! .. 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, ifail, lda, ldb, lenc, lmax, &
locc, m, n, nf, nwcm, nwcn, nwct, &
nwl, want_coeffs, want_level
Character (10) :: mode, wavnam, wtrans
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: a(:,:), b(:,:), c(:), d(:,:)
Integer, Allocatable :: dwtlevm(:), dwtlevn(:)
Integer :: icomm(180)
! .. Intrinsic Procedures ..
Intrinsic :: sum
! .. Executable Statements ..
Continue
Write (nout,*) 'C09ABF Example Program Results'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
! Read problem parameters
Read (nin,*) m, n
lda = m
ldb = m
Read (nin,*) wavnam, mode
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: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call c09abf(wavnam,wtrans,mode,m,n,lmax,nf,nwct,nwcn,icomm,ifail)
lenc = nwct
Allocate (c(lenc),dwtlevm(lmax),dwtlevn(lmax))
! Calculate one less than the max possible number of levels
nwl = lmax - 1
! Perform Discrete Wavelet transform
ifail = 0
Call c09ecf(m,n,a,lda,lenc,c,nwl,dwtlevm,dwtlevn,icomm,ifail)
! c09abf returns nwct based on max levels, so recalculate.
nwct = sum(3*dwtlevm(1:nwl)*dwtlevn(1:nwl))
nwct = nwct + dwtlevm(1)*dwtlevn(1)
! Print the details of the transform.
Write (nout,*)
Write (nout,99997) nf
Write (nout,99996) nwl
Write (nout,99995)
Write (nout,99994) dwtlevm(1:nwl)
Write (nout,99993)
Write (nout,99994) dwtlevn(1:nwl)
Write (nout,99992) nwct
Write (nout,*)
Write (nout,99991)
Write (nout,99998) c(1:nwct)
! Now select a nominated matrix of coefficients at a nominated level.
! Remember that level 0 is input data, 1 first coeffs and so on up to nwl, which
! is the deepest level and contains approx. coefficients.
want_level = nwl - 1
! Print only veritcal detail coeffs at selected level.
want_coeffs = 1
nwcm = dwtlevm(nwl-want_level+1)
nwcn = dwtlevn(nwl-want_level+1)
Allocate (d(nwcm,nwcn))
! Extract the selected set of coefficients. Firstly figure out
! where in C we need to copy from.
!
! Count number of coefficients at level 1 to skip past. Remember level
! 1 contains 1 set of approximation and 3 sets of detail coefficients.
locc = 4*dwtlevm(1)*dwtlevn(1)
! Count number of coefficients at levels 2 up to want_level-1, which
! each contain 3 sets of detail coefficients.
Do i = 2, want_level - 1
locc = 3*dwtlevm(i)*dwtlevn(i)
End Do
! Now locc points to the level we want.
! Coefficients are stored in the order: vertical,
! horizontal, diagonal. We want vertical, which is first:
locc = locc + 1
! Extract the coefficients from c starting at locc
Do i = 1, nwcn
d(1:nwcm,i) = c(locc:locc+nwcm-1)
locc = locc + nwcm
End Do
! Print out the selected coefficients
Write (nout,*)
Write (nout,99989) want_coeffs, want_level
Do i = 1, nwcm
Write (nout,99998) d(i,1:nwcn)
End Do
! Reconstruct original data
ifail = 0
Call c09edf(nwl,lenc,c,m,n,b,ldb,icomm,ifail)
Write (nout,*)
Write (nout,99990)
Do i = 1, m
Write (nout,99998) b(i,1:n)
End Do
99999 Format (1X,' Parameters read from file :: '/' Wavelet : ',A10, &
' End mode : ',A10,' M = ',I10,' N = ',I10)
99998 Format (8(F8.4,1X):)
99997 Format (1X,' Length of wavelet filter : ',I10)
99996 Format (1X,' Number of Levels : ',I10)
99995 Format (1X, &
' Number of coefficients in first dimension for each level : ')
99994 Format (16X,8(I8,1X):)
99993 Format (1X, &
' Number of coefficients in second dimension for each level : ')
99992 Format (1X,' Total number of wavelet coefficients : ',I10)
99991 Format (1X,' Wavelet coefficients C : ')
99990 Format (1X,' Reconstruction B : ')
99989 Format (1X,' Type ',I1,' coefficients at selected wavelet level, ',I4, &
': ')
End Program c09abfe