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:   PetscRandom rdm;           /* for creating random input */
 19:   PetscScalar a;             /* used to scale output */
 20:   PetscReal   enorm;         /* norm for sanity check */
 21:   PetscInt    n = 10, N = 1; /* FFT dimension params */
 22:   PetscInt    DIM, dim[5];   /* FFT params */

 24:   PetscFunctionBeginUser;
 25:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 26:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));

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

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

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

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

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

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

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

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

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

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

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

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

 91: /*TEST

 93:     build:
 94:       requires: fftw complex

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

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

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

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

116: TEST*/