1 static char help[] = "Test sequential FFTW convolution\n\n"; 2 3 /* 4 Compiling the code: 5 This code uses the complex numbers, so configure must be given --with-scalar-type=complex to enable this 6 */ 7 8 #include <petscmat.h> 9 10 int main(int argc, char **args) 11 { 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, j; 22 Vec w, x, y1, y2, z1, z2; 23 PetscScalar *a, *a2, *a3; 24 PetscScalar s; 25 PetscRandom rdm; 26 PetscReal enorm; 27 PetscInt func = 0; 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", "ex121", 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 /* create vectors of length N=n^DIM */ 47 for (i = 0, N = 1; i < DIM; i++) N *= dim[i]; 48 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n", DIM, N)); 49 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x)); 50 PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector")); 51 PetscCall(VecDuplicate(x, &w)); 52 PetscCall(PetscObjectSetName((PetscObject)w, "Window vector")); 53 PetscCall(VecDuplicate(x, &y1)); 54 PetscCall(PetscObjectSetName((PetscObject)y1, "Frequency space vector")); 55 PetscCall(VecDuplicate(x, &y2)); 56 PetscCall(PetscObjectSetName((PetscObject)y2, "Frequency space window vector")); 57 PetscCall(VecDuplicate(x, &z1)); 58 PetscCall(PetscObjectSetName((PetscObject)z1, "Reconstructed convolution")); 59 PetscCall(VecDuplicate(x, &z2)); 60 PetscCall(PetscObjectSetName((PetscObject)z2, "Real space convolution")); 61 62 if (function == RANDOM) { 63 PetscCall(VecSetRandom(x, rdm)); 64 } else if (function == CONSTANT) { 65 PetscCall(VecSet(x, 1.0)); 66 } else if (function == TANH) { 67 PetscCall(VecGetArray(x, &a)); 68 for (i = 0; i < N; ++i) a[i] = tanh((i - N / 2.0) * (10.0 / N)); 69 PetscCall(VecRestoreArray(x, &a)); 70 } 71 if (view) PetscCall(VecView(x, PETSC_VIEWER_DRAW_WORLD)); 72 73 /* Create window function */ 74 PetscCall(VecGetArray(w, &a)); 75 for (i = 0; i < N; ++i) { 76 /* Step Function */ 77 a[i] = (i > N / 4 && i < 3 * N / 4) ? 1.0 : 0.0; 78 /* Delta Function */ 79 /*a[i] = (i == N/2)? 1.0: 0.0; */ 80 } 81 PetscCall(VecRestoreArray(w, &a)); 82 if (view) PetscCall(VecView(w, PETSC_VIEWER_DRAW_WORLD)); 83 84 /* create FFTW object */ 85 PetscCall(MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A)); 86 87 /* Convolve x with w*/ 88 PetscCall(MatMult(A, x, y1)); 89 PetscCall(MatMult(A, w, y2)); 90 PetscCall(VecPointwiseMult(y1, y1, y2)); 91 if (view && i == 0) PetscCall(VecView(y1, PETSC_VIEWER_DRAW_WORLD)); 92 PetscCall(MatMultTranspose(A, y1, z1)); 93 94 /* Compute the real space convolution */ 95 PetscCall(VecGetArray(x, &a)); 96 PetscCall(VecGetArray(w, &a2)); 97 PetscCall(VecGetArray(z2, &a3)); 98 for (i = 0; i < N; ++i) { 99 /* PetscInt checkInd = (i > N/2-1)? i-N/2: i+N/2;*/ 100 101 a3[i] = 0.0; 102 for (j = -N / 2 + 1; j < N / 2; ++j) { 103 PetscInt xpInd = (j < 0) ? N + j : j; 104 PetscInt diffInd = (i - j < 0) ? N - (j - i) : (i - j > N - 1) ? i - j - N : i - j; 105 106 a3[i] += a[xpInd] * a2[diffInd]; 107 } 108 } 109 PetscCall(VecRestoreArray(x, &a)); 110 PetscCall(VecRestoreArray(w, &a2)); 111 PetscCall(VecRestoreArray(z2, &a3)); 112 113 /* compare z1 and z2. FFTW computes an unnormalized DFT, thus z1 = N*z2 */ 114 s = 1.0 / (PetscReal)N; 115 PetscCall(VecScale(z1, s)); 116 if (view) PetscCall(VecView(z1, PETSC_VIEWER_DRAW_WORLD)); 117 if (view) PetscCall(VecView(z2, PETSC_VIEWER_DRAW_WORLD)); 118 PetscCall(VecAXPY(z1, -1.0, z2)); 119 PetscCall(VecNorm(z1, NORM_1, &enorm)); 120 if (enorm > 1.e-11) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |z1 - z2| %g\n", (double)enorm)); 121 122 /* free spaces */ 123 PetscCall(VecDestroy(&x)); 124 PetscCall(VecDestroy(&y1)); 125 PetscCall(VecDestroy(&y2)); 126 PetscCall(VecDestroy(&z1)); 127 PetscCall(VecDestroy(&z2)); 128 PetscCall(VecDestroy(&w)); 129 PetscCall(MatDestroy(&A)); 130 } 131 PetscCall(PetscRandomDestroy(&rdm)); 132 PetscCall(PetscFinalize()); 133 return 0; 134 } 135 136 /*TEST 137 138 build: 139 requires: fftw complex 140 141 test: 142 output_file: output/ex121.out 143 TODO: Example or FFTW interface is broken 144 145 TEST*/ 146