/* nag_mldwt_3d (c09fcc) Example Program. * * Copyright 2013 Numerical Algorithms Group. * * Mark 24, 2013. */ #include #include #include #include #include #include #define A(I,J,K) a[I-1 + (J-1)* lda + (K-1)* lda * sda] #define B(I,J,K) b[I-1 + (J-1)* ldb + (K-1)* ldb * sdb] #define E(I,J,K) e[I-1 + (J-1)* m + (K-1)* m * n] #define D(I,J,K) d[I-1 + (J-1)* ldd + (K-1)* ldd * sdd] int main(void) { /* Scalars */ Integer exit_status = 0; Integer lda, ldb, ldd, sda, sdb, sdd, lenc, i, j, k; Integer m, n, fr, nwcfr, nwcm, nwcn, nwct, nwlmax, nwl, nwlinv, nf; Integer want_coeffs, want_level; double eps, esq, frob; /* Arrays */ char mode[25], wavnam[25]; double *a = 0, *b = 0, *c = 0, *d = 0, *e = 0; Integer *dwtlvfr = 0, *dwtlvm = 0, *dwtlvn = 0; Integer icomm[260]; /* Nag Types */ Nag_Wavelet wavnamenum; Nag_WaveletMode modenum; Nag_MatrixType matrix = Nag_GeneralMatrix; Nag_OrderType order = Nag_ColMajor; Nag_DiagType diag = Nag_NonUnitDiag; NagError fail; INIT_FAIL(fail); printf("nag_mldwt_3d (c09fcc) Example Program Results\n\n"); /* Skip heading in data file and read problem parameters */ scanf("%*[^\n] %ld%ld%ld%*[^\n]", &m, &n, &fr); lda = m; ldb = m; sda = n; sdb = n; scanf("%25s%25s%*[^\n]\n", wavnam, mode); if (!(a = NAG_ALLOC((lda)*(sda)*(fr), double)) || !(b = NAG_ALLOC((ldb)*(sdb)*(fr), double))|| !(e = NAG_ALLOC((m)*(n)*(fr), double))) { printf("Allocation failure\n"); exit_status = 1; goto END; } printf("Parameters read from file :: \n"); printf("MLDWT :: Wavelet : %s\n", wavnam); printf(" End mode : %s\n", mode); printf(" m : %4ld\n", m); printf(" n : %4ld\n", n); printf(" fr : %4ld\n\n", fr); /* nag_enum_name_to_value (x04nac). * Converts NAG enum member name to value */ wavnamenum = (Nag_Wavelet) nag_enum_name_to_value(wavnam); modenum = (Nag_WaveletMode) nag_enum_name_to_value(mode); /* Read data array */ for (k=1; k<=fr; k++) { for (i=1; i<=m; i++) { for (j=1; j<=n; j++) scanf("%lf", &A(i, j, k)); scanf("%*[^\n] "); } scanf("%*[^\n] "); } /* Print out the input data */ printf("Input Data :\n"); for (k=1; k<=fr; k++) { /* nag_gen_real_mat_print_comp (x04cbc). * Prints out a matrix. */ nag_gen_real_mat_print_comp(order, matrix, diag, m, n, &A(1, 1, k), lda, "%8.4f"," ", Nag_NoLabels, 0, Nag_NoLabels, 0, 80, 0, 0, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n", fail.message); exit_status = 2; goto END; } printf("\n"); } /* nag_wfilt_3d (c09acc). * Three-dimensional wavelet filter initialization */ nag_wfilt_3d(wavnamenum, Nag_MultiLevel, modenum, m, n, fr, &nwlmax, &nf, &nwct, &nwcn, &nwcfr, icomm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_wfilt_3d (c09acc).\n%s\n", fail.message); exit_status = 3; goto END; } lenc = nwct; if (!(c = NAG_ALLOC((lenc), double)) || !(dwtlvm = NAG_ALLOC((nwlmax), Integer)) || !(dwtlvn = NAG_ALLOC((nwlmax), Integer)) || !(dwtlvfr = NAG_ALLOC((nwlmax), Integer))) { printf("Allocation failure\n"); exit_status = 4; goto END; } nwl = nwlmax; /* nag_mldwt_3d (c09fcc). * Three-dimensional multi-level discrete wavelet transform */ nag_mldwt_3d(m, n, fr, a, lda, sda, lenc, c, nwl, dwtlvm, dwtlvn, dwtlvfr, icomm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_mldwt_3d (c09fcc).\n%s\n", fail.message); exit_status = 5; goto END; } printf("Number of Levels : %4ld\n", nwl); printf("Number of coefficients in 1st dimension for each level:\n"); for (i=0; ieps) { printf("\nFail: Frobenius norm of B-A, where A is the original \n" "data and B is the reconstrucion, is too large.\n"); } else { printf("\nSuccess: the reconstruction matches the original.\n"); } END: NAG_FREE(a); NAG_FREE(b); NAG_FREE(c); NAG_FREE(d); NAG_FREE(e); NAG_FREE(dwtlvfr); NAG_FREE(dwtlvm); NAG_FREE(dwtlvn); return exit_status; }