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

 29:   PetscInitialize(&argc,&args,(char*)0,help);
 30: #if !defined(PETSC_USE_COMPLEX)
 31:   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires complex numbers");
 32: #endif
 33:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 34:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This is a uniprocessor example only!");
 35:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex112");
 36:   PetscOptionsEList("-function", "Function type", "ex112", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
 37:   PetscOptionsBool("-vec_view_draw", "View the functions", "ex112", view, &view, NULL);
 38:   function = (FuncType) func;
 39:   PetscOptionsEnd();

 41:   for (DIM = 0; DIM < ndim; DIM++) {
 42:     dim[DIM] = n;  /* size of transformation in DIM-dimension */
 43:   }
 44:   PetscRandomCreate(PETSC_COMM_SELF, &rdm);
 45:   PetscRandomSetFromOptions(rdm);

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

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

 54:     /* create vectors of length N=n^DIM */
 55:     MatGetVecs(A,&x,&y);
 56:     MatGetVecs(A,&z,NULL);
 57:     PetscObjectSetName((PetscObject) x, "Real space vector");
 58:     PetscObjectSetName((PetscObject) y, "Frequency space vector");
 59:     PetscObjectSetName((PetscObject) z, "Reconstructed vector");

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

 76:     /* apply FFTW_FORWARD and FFTW_BACKWARD several times on same x, y, and z */
 77:     for (i=0; i<3; i++) {
 78:       MatMult(A,x,y);
 79:       if (view && i == 0) {VecView(y, PETSC_VIEWER_DRAW_WORLD);}
 80:       MatMultTranspose(A,y,z);

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

 93:     /* apply FFTW_FORWARD and FFTW_BACKWARD several times on different x */
 94:     for (i=0; i<3; i++) {
 95:       VecDestroy(&x);
 96:       VecCreateSeq(PETSC_COMM_SELF,N,&x);
 97:       VecSetRandom(x, rdm);

 99:       MatMult(A,x,y);
100:       MatMultTranspose(A,y,z);

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

113:     /* free spaces */
114:     VecDestroy(&x);
115:     VecDestroy(&y);
116:     VecDestroy(&z);
117:     MatDestroy(&A);
118:   }
119:   PetscRandomDestroy(&rdm);
120:   PetscFinalize();
121:   return 0;
122: }