Actual source code: ex228.c

  1: static char help[] = "Test duplication/destruction of FFTW vecs \n\n";

  3: /*
  4:  Compiling the code:
  5:    This code uses the FFTW interface.
  6:    Use one of the options below to configure:
  7:    --with-fftw-dir=/.... or --download-fftw
  8:  Usage:
  9:    mpiexec -np <np> ./ex228
 10: */

 12: #include <petscmat.h>
 13: int main(int argc, char **args)
 14: {
 15:   Mat         A;             /* FFT Matrix */
 16:   Vec         x, y, z;       /* Work vectors */
 17:   Vec         x1, y1, z1;    /* Duplicate vectors */
 18:   PetscInt    i, k;          /* for iterating over dimensions */
 19:   PetscRandom rdm;           /* for creating random input */
 20:   PetscScalar a;             /* used to scale output */
 21:   PetscReal   enorm;         /* norm for sanity check */
 22:   PetscInt    n = 10, N = 1; /* FFT dimension params */
 23:   PetscInt    DIM, dim[5];   /* FFT params */

 26:   PetscInitialize(&argc, &args, (char *)0, help);
 27:   PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);

 29:   /* To create random input vector */
 30:   PetscRandomCreate(PETSC_COMM_SELF, &rdm);
 31:   PetscRandomSetFromOptions(rdm);

 33:   /* Iterate over dimensions, use PETSc-FFTW interface */
 34:   for (i = 1; i < 5; i++) {
 35:     DIM = i;
 36:     N   = 1;
 37:     for (k = 0; k < i; k++) {
 38:       dim[k] = n;
 39:       N *= n;
 40:     }

 42:     PetscPrintf(PETSC_COMM_WORLD, "\n %" PetscInt_FMT " dimensions: FFTW on vector of size %" PetscInt_FMT " \n", DIM, N);

 44:     /* create FFTW object */
 45:     MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A);
 46:     /* create vectors of length N */
 47:     MatCreateVecsFFTW(A, &x, &y, &z);

 49:     PetscObjectSetName((PetscObject)x, "Real space vector");
 50:     PetscObjectSetName((PetscObject)y, "Frequency space vector");
 51:     PetscObjectSetName((PetscObject)z, "Reconstructed vector");

 53:     /* Test vector duplication*/
 54:     VecDuplicate(x, &x1);
 55:     VecDuplicate(y, &y1);
 56:     VecDuplicate(z, &z1);

 58:     /* Set values of space vector x, copy to duplicate */
 59:     VecSetRandom(x, rdm);
 60:     VecCopy(x, x1);

 62:     /* Apply FFTW_FORWARD and FFTW_BACKWARD */
 63:     MatMult(A, x, y);
 64:     MatMultTranspose(A, y, z);

 66:     /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */
 67:     MatMult(A, x1, y1);
 68:     MatMultTranspose(A, y1, z1);

 70:     /* Compare x and z1. FFTW computes an unnormalized DFT, thus z1 = N*x */
 71:     a = 1.0 / (PetscReal)N;
 72:     VecScale(z1, a);
 73:     VecAXPY(z1, -1.0, x);
 74:     VecNorm(z1, NORM_1, &enorm);
 75:     if (enorm > 1.e-9) PetscPrintf(PETSC_COMM_WORLD, "  Error norm of |x - z1| %g\n", enorm);

 77:     /* free spaces */
 78:     VecDestroy(&x1);
 79:     VecDestroy(&y1);
 80:     VecDestroy(&z1);
 81:     VecDestroy(&x);
 82:     VecDestroy(&y);
 83:     VecDestroy(&z);
 84:     MatDestroy(&A);
 85:   }

 87:   PetscRandomDestroy(&rdm);
 88:   PetscFinalize();
 89:   return 0;
 90: }

 92: /*TEST

 94:     build:
 95:       requires: fftw complex

 97:     test:
 98:       suffix: 2
 99:       nsize : 4
100:       args: -mat_fftw_plannerflags FFTW_ESTIMATE -n 16

102:     test:
103:       suffix: 3
104:       nsize : 2
105:       args: -mat_fftw_plannerflags FFTW_MEASURE -n 12

107:     test:
108:       suffix: 4
109:       nsize : 2
110:       args: -mat_fftw_plannerflags FFTW_PATIENT -n 10

112:     test:
113:       suffix: 5
114:       nsize : 1
115:       args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE -n 5

117: TEST*/