Program c09ecfe ! 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, 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) ! Extract coefficients 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