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

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

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

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

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

 40:     /* create FFTW object */
 41:     MatCreateFFT(PETSC_COMM_SELF,DIM,dim,MATFFTW,&A);
 42:     /* create vectors of length N */
 43:     MatCreateVecsFFTW(A,&x,&y,&z);

 45:     PetscObjectSetName((PetscObject) x, "Real space vector");
 46:     PetscObjectSetName((PetscObject) y, "Frequency space vector");
 47:     PetscObjectSetName((PetscObject) z, "Reconstructed vector");

 49:     /* Test vector duplication*/
 50:     VecDuplicate(x,&x1);
 51:     VecDuplicate(y,&y1);
 52:     VecDuplicate(z,&z1);

 54:     /* Set values of space vector x, copy to duplicate */
 55:     VecSetRandom(x,rdm);
 56:     VecCopy(x,x1);

 58:     /* Apply FFTW_FORWARD and FFTW_BACKWARD */
 59:     MatMult(A,x,y);
 60:     MatMultTranspose(A,y,z);

 62:     /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */
 63:     MatMult(A,x1,y1);
 64:     MatMultTranspose(A,y1,z1);

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

 73:     /* free spaces */
 74:     VecDestroy(&x1);
 75:     VecDestroy(&y1);
 76:     VecDestroy(&z1);
 77:     VecDestroy(&x);
 78:     VecDestroy(&y);
 79:     VecDestroy(&z);
 80:     MatDestroy(&A);
 81:   }

 83:   PetscRandomDestroy(&rdm);
 84:   PetscFinalize();
 85:   return 0;
 86: }

 88: /*TEST

 90:     build:
 91:       requires: fftw complex

 93:     test:
 94:       suffix: 2
 95:       nsize : 4
 96:       args: -mat_fftw_plannerflags FFTW_ESTIMATE -n 16

 98:     test:
 99:       suffix: 3
100:       nsize : 2
101:       args: -mat_fftw_plannerflags FFTW_MEASURE -n 12

103:     test:
104:       suffix: 4
105:       nsize : 2
106:       args: -mat_fftw_plannerflags FFTW_PATIENT -n 10

108:     test:
109:       suffix: 5
110:       nsize : 1
111:       args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE -n 5

113: TEST*/