Actual source code: ex112.c

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

  3: /*
  4:   Compiling the code:
  5:       This code uses the complex numbers version of PETSc, so configure
  6:       must be run to enable this

  8: */

 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;

 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", "ex112", 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:     for (i = 0, N = 1; i < DIM; i++) N *= dim[i];
 47:     PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n", DIM, N);

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

 52:     /* create vectors of length N=n^DIM */
 53:     MatCreateVecs(A, &x, &y);
 54:     MatCreateVecs(A, &z, NULL);
 55:     PetscObjectSetName((PetscObject)x, "Real space vector");
 56:     PetscObjectSetName((PetscObject)y, "Frequency space vector");
 57:     PetscObjectSetName((PetscObject)z, "Reconstructed vector");

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

 72:     /* apply FFTW_FORWARD and FFTW_BACKWARD several times on same x, y, and z */
 73:     for (i = 0; i < 3; i++) {
 74:       MatMult(A, x, y);
 75:       if (view && i == 0) VecView(y, PETSC_VIEWER_DRAW_WORLD);
 76:       MatMultTranspose(A, y, z);

 78:       /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
 79:       s = 1.0 / (PetscReal)N;
 80:       VecScale(z, s);
 81:       if (view && i == 0) VecView(z, PETSC_VIEWER_DRAW_WORLD);
 82:       VecAXPY(z, -1.0, x);
 83:       VecNorm(z, NORM_1, &enorm);
 84:       if (enorm > tol) PetscPrintf(PETSC_COMM_SELF, "  Error norm of |x - z| %g\n", (double)enorm);
 85:     }

 87:     /* apply FFTW_FORWARD and FFTW_BACKWARD several times on different x */
 88:     for (i = 0; i < 3; i++) {
 89:       VecDestroy(&x);
 90:       VecCreateSeq(PETSC_COMM_SELF, N, &x);
 91:       VecSetRandom(x, rdm);

 93:       MatMult(A, x, y);
 94:       MatMultTranspose(A, y, z);

 96:       /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
 97:       s = 1.0 / (PetscReal)N;
 98:       VecScale(z, s);
 99:       if (view && i == 0) VecView(z, PETSC_VIEWER_DRAW_WORLD);
100:       VecAXPY(z, -1.0, x);
101:       VecNorm(z, NORM_1, &enorm);
102:       if (enorm > tol) PetscPrintf(PETSC_COMM_SELF, "  Error norm of new |x - z| %g\n", (double)enorm);
103:     }

105:     /* free spaces */
106:     VecDestroy(&x);
107:     VecDestroy(&y);
108:     VecDestroy(&z);
109:     MatDestroy(&A);
110:   }
111:   PetscRandomDestroy(&rdm);
112:   PetscFinalize();
113:   return 0;
114: }

116: /*TEST

118:    build:
119:       requires:  fftw complex

121:    test:
122:       args: -mat_fftw_plannerflags FFTW_ESTIMATE
123:       output_file: output/ex112.out

125:    test:
126:       suffix: 2
127:       args: -mat_fftw_plannerflags FFTW_MEASURE
128:       output_file: output/ex112.out
129:       requires: !defined(PETSC_USE_CXXCOMPLEX)

131:    test:
132:       suffix: 3
133:       args: -mat_fftw_plannerflags FFTW_PATIENT
134:       output_file: output/ex112.out

136:    test:
137:       suffix: 4
138:       args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE
139:       output_file: output/ex112.out

141: TEST*/