Actual source code: ex121.c

  1: static char help[] = "Test sequential FFTW convolution\n\n";

  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: */

  8: #include <petscmat.h>

 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;

 32:   PetscInitialize(&argc, &args, (char *)0, help);
 33:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 35:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex112");
 36:   PetscOptionsEList("-function", "Function type", "ex121", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
 37:   PetscOptionsBool("-vec_view draw", "View the functions", "ex112", view, &view, NULL);
 38:   function = (FuncType)func;
 39:   PetscOptionsEnd();

 41:   for (DIM = 0; DIM < ndim; DIM++) { dim[DIM] = n; /* size of transformation in DIM-dimension */ }
 42:   PetscRandomCreate(PETSC_COMM_SELF, &rdm);
 43:   PetscRandomSetFromOptions(rdm);

 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:     PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n", DIM, N);
 49:     VecCreateSeq(PETSC_COMM_SELF, N, &x);
 50:     PetscObjectSetName((PetscObject)x, "Real space vector");
 51:     VecDuplicate(x, &w);
 52:     PetscObjectSetName((PetscObject)w, "Window vector");
 53:     VecDuplicate(x, &y1);
 54:     PetscObjectSetName((PetscObject)y1, "Frequency space vector");
 55:     VecDuplicate(x, &y2);
 56:     PetscObjectSetName((PetscObject)y2, "Frequency space window vector");
 57:     VecDuplicate(x, &z1);
 58:     PetscObjectSetName((PetscObject)z1, "Reconstructed convolution");
 59:     VecDuplicate(x, &z2);
 60:     PetscObjectSetName((PetscObject)z2, "Real space convolution");

 62:     if (function == RANDOM) {
 63:       VecSetRandom(x, rdm);
 64:     } else if (function == CONSTANT) {
 65:       VecSet(x, 1.0);
 66:     } else if (function == TANH) {
 67:       VecGetArray(x, &a);
 68:       for (i = 0; i < N; ++i) a[i] = tanh((i - N / 2.0) * (10.0 / N));
 69:       VecRestoreArray(x, &a);
 70:     }
 71:     if (view) VecView(x, PETSC_VIEWER_DRAW_WORLD);

 73:     /* Create window function */
 74:     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:     VecRestoreArray(w, &a);
 82:     if (view) VecView(w, PETSC_VIEWER_DRAW_WORLD);

 84:     /* create FFTW object */
 85:     MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A);

 87:     /* Convolve x with w*/
 88:     MatMult(A, x, y1);
 89:     MatMult(A, w, y2);
 90:     VecPointwiseMult(y1, y1, y2);
 91:     if (view && i == 0) VecView(y1, PETSC_VIEWER_DRAW_WORLD);
 92:     MatMultTranspose(A, y1, z1);

 94:     /* Compute the real space convolution */
 95:     VecGetArray(x, &a);
 96:     VecGetArray(w, &a2);
 97:     VecGetArray(z2, &a3);
 98:     for (i = 0; i < N; ++i) {
 99:       /* PetscInt checkInd = (i > N/2-1)? i-N/2: i+N/2;*/

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;

106:         a3[i] += a[xpInd] * a2[diffInd];
107:       }
108:     }
109:     VecRestoreArray(x, &a);
110:     VecRestoreArray(w, &a2);
111:     VecRestoreArray(z2, &a3);

113:     /* compare z1 and z2. FFTW computes an unnormalized DFT, thus z1 = N*z2 */
114:     s = 1.0 / (PetscReal)N;
115:     VecScale(z1, s);
116:     if (view) VecView(z1, PETSC_VIEWER_DRAW_WORLD);
117:     if (view) VecView(z2, PETSC_VIEWER_DRAW_WORLD);
118:     VecAXPY(z1, -1.0, z2);
119:     VecNorm(z1, NORM_1, &enorm);
120:     if (enorm > 1.e-11) PetscPrintf(PETSC_COMM_SELF, "  Error norm of |z1 - z2| %g\n", (double)enorm);

122:     /* free spaces */
123:     VecDestroy(&x);
124:     VecDestroy(&y1);
125:     VecDestroy(&y2);
126:     VecDestroy(&z1);
127:     VecDestroy(&z2);
128:     VecDestroy(&w);
129:     MatDestroy(&A);
130:   }
131:   PetscRandomDestroy(&rdm);
132:   PetscFinalize();
133:   return 0;
134: }

136: /*TEST

138:    build:
139:       requires: fftw complex

141:    test:
142:       output_file: output/ex121.out
143:       TODO: Example or FFTW interface is broken

145: TEST*/