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 {RANDOM, CONSTANT, TANH, NUM_FUNCS} FuncType;
 14:   const char     *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
 15:   Mat            A;
 16:   PetscMPIInt    size;
 17:   PetscInt       n = 10,N,ndim=4,dim[4],DIM,i;
 18:   Vec            x,y,z;
 19:   PetscScalar    s;
 20:   PetscRandom    rdm;
 21:   PetscReal      enorm, tol = PETSC_SMALL;
 22:   PetscInt       func;
 23:   FuncType       function = RANDOM;
 24:   PetscBool      view     = PETSC_FALSE;

 27:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 28:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 29:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This is a uniprocessor example only!");
 30:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex112");
 31:   PetscOptionsEList("-function", "Function type", "ex112", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
 32:   PetscOptionsBool("-vec_view_draw", "View the functions", "ex112", view, &view, NULL);
 33:   function = (FuncType) func;
 34:   PetscOptionsEnd();

 36:   for (DIM = 0; DIM < ndim; DIM++) {
 37:     dim[DIM] = n;  /* size of transformation in DIM-dimension */
 38:   }
 39:   PetscRandomCreate(PETSC_COMM_SELF, &rdm);
 40:   PetscRandomSetFromOptions(rdm);

 42:   for (DIM = 1; DIM < 5; DIM++) {
 43:     for (i = 0, N = 1; i < DIM; i++) N *= dim[i];
 44:     PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n",DIM,N);

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

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

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

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

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

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

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

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

108:     /* free spaces */
109:     VecDestroy(&x);
110:     VecDestroy(&y);
111:     VecDestroy(&z);
112:     MatDestroy(&A);
113:   }
114:   PetscRandomDestroy(&rdm);
115:   PetscFinalize();
116:   return ierr;
117: }


120: /*TEST

122:    build:
123:       requires:  fftw complex

125:    test:
126:       args: -mat_fftw_plannerflags FFTW_ESTIMATE
127:       output_file: output/ex112.out

129:    test:
130:       suffix: 2
131:       args: -mat_fftw_plannerflags FFTW_MEASURE
132:       output_file: output/ex112.out
133:       requires: !define(PETSC_USE_CXXCOMPLEX)

135:    test:
136:       suffix: 3
137:       args: -mat_fftw_plannerflags FFTW_PATIENT
138:       output_file: output/ex112.out

140:    test:
141:       suffix: 4
142:       args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE
143:       output_file: output/ex112.out

145: TEST*/