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;

 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();

 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));

 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));

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

 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"));

 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));

 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));

 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:     }

 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));

 93:       PetscCall(MatMult(A, x, y));
 94:       PetscCall(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:       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:     }

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: }

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