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