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