1 static char help[] = "Test sequential FFTW interface \n\n"; 2 3 /* 4 Compiling the code: 5 This code uses the complex numbers version of PETSc, so configure 6 must be run to enable this 7 8 */ 9 10 #include <petscmat.h> 11 int main(int argc, char **args) { 12 typedef enum { 13 RANDOM, 14 CONSTANT, 15 TANH, 16 NUM_FUNCS 17 } FuncType; 18 const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"}; 19 Mat A; 20 PetscMPIInt size; 21 PetscInt n = 10, N, ndim = 4, dim[4], DIM, i; 22 Vec x, y, z; 23 PetscScalar s; 24 PetscRandom rdm; 25 PetscReal enorm, tol = PETSC_SMALL; 26 PetscInt func; 27 FuncType function = RANDOM; 28 PetscBool view = PETSC_FALSE; 29 30 PetscFunctionBeginUser; 31 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 32 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 33 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 34 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex112"); 35 PetscCall(PetscOptionsEList("-function", "Function type", "ex112", funcNames, NUM_FUNCS, funcNames[function], &func, NULL)); 36 PetscCall(PetscOptionsBool("-vec_view_draw", "View the functions", "ex112", view, &view, NULL)); 37 function = (FuncType)func; 38 PetscOptionsEnd(); 39 40 for (DIM = 0; DIM < ndim; DIM++) { dim[DIM] = n; /* size of transformation in DIM-dimension */ } 41 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm)); 42 PetscCall(PetscRandomSetFromOptions(rdm)); 43 44 for (DIM = 1; DIM < 5; DIM++) { 45 for (i = 0, N = 1; i < DIM; i++) N *= dim[i]; 46 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n", DIM, N)); 47 48 /* create FFTW object */ 49 PetscCall(MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A)); 50 51 /* create vectors of length N=n^DIM */ 52 PetscCall(MatCreateVecs(A, &x, &y)); 53 PetscCall(MatCreateVecs(A, &z, NULL)); 54 PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector")); 55 PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector")); 56 PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector")); 57 58 /* set values of space vector x */ 59 if (function == RANDOM) { 60 PetscCall(VecSetRandom(x, rdm)); 61 } else if (function == CONSTANT) { 62 PetscCall(VecSet(x, 1.0)); 63 } else if (function == TANH) { 64 PetscScalar *a; 65 PetscCall(VecGetArray(x, &a)); 66 for (i = 0; i < N; ++i) a[i] = tanh((i - N / 2.0) * (10.0 / N)); 67 PetscCall(VecRestoreArray(x, &a)); 68 } 69 if (view) PetscCall(VecView(x, PETSC_VIEWER_DRAW_WORLD)); 70 71 /* apply FFTW_FORWARD and FFTW_BACKWARD several times on same x, y, and z */ 72 for (i = 0; i < 3; i++) { 73 PetscCall(MatMult(A, x, y)); 74 if (view && i == 0) PetscCall(VecView(y, PETSC_VIEWER_DRAW_WORLD)); 75 PetscCall(MatMultTranspose(A, y, z)); 76 77 /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */ 78 s = 1.0 / (PetscReal)N; 79 PetscCall(VecScale(z, s)); 80 if (view && i == 0) PetscCall(VecView(z, PETSC_VIEWER_DRAW_WORLD)); 81 PetscCall(VecAXPY(z, -1.0, x)); 82 PetscCall(VecNorm(z, NORM_1, &enorm)); 83 if (enorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %g\n", (double)enorm)); 84 } 85 86 /* apply FFTW_FORWARD and FFTW_BACKWARD several times on different x */ 87 for (i = 0; i < 3; i++) { 88 PetscCall(VecDestroy(&x)); 89 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x)); 90 PetscCall(VecSetRandom(x, rdm)); 91 92 PetscCall(MatMult(A, x, y)); 93 PetscCall(MatMultTranspose(A, y, z)); 94 95 /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */ 96 s = 1.0 / (PetscReal)N; 97 PetscCall(VecScale(z, s)); 98 if (view && i == 0) PetscCall(VecView(z, PETSC_VIEWER_DRAW_WORLD)); 99 PetscCall(VecAXPY(z, -1.0, x)); 100 PetscCall(VecNorm(z, NORM_1, &enorm)); 101 if (enorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of new |x - z| %g\n", (double)enorm)); 102 } 103 104 /* free spaces */ 105 PetscCall(VecDestroy(&x)); 106 PetscCall(VecDestroy(&y)); 107 PetscCall(VecDestroy(&z)); 108 PetscCall(MatDestroy(&A)); 109 } 110 PetscCall(PetscRandomDestroy(&rdm)); 111 PetscCall(PetscFinalize()); 112 return 0; 113 } 114 115 /*TEST 116 117 build: 118 requires: fftw complex 119 120 test: 121 args: -mat_fftw_plannerflags FFTW_ESTIMATE 122 output_file: output/ex112.out 123 124 test: 125 suffix: 2 126 args: -mat_fftw_plannerflags FFTW_MEASURE 127 output_file: output/ex112.out 128 requires: !defined(PETSC_USE_CXXCOMPLEX) 129 130 test: 131 suffix: 3 132 args: -mat_fftw_plannerflags FFTW_PATIENT 133 output_file: output/ex112.out 134 135 test: 136 suffix: 4 137 args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE 138 output_file: output/ex112.out 139 140 TEST*/ 141