Actual source code: ex112.c

  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: PetscInt main(PetscInt 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:   PetscTruth     view = PETSC_FALSE;

 29:   PetscInitialize(&argc,&args,(char *)0,help);
 30: #if !defined(PETSC_USE_COMPLEX)
 31:   SETERRQ(PETSC_ERR_SUP, "This example requires complex numbers");
 32: #endif
 33:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 34:   if (size != 1) SETERRQ(PETSC_ERR_SUP, "This is a uniprocessor example only!");
 35:   PetscOptionsBegin(PETSC_COMM_WORLD, PETSC_NULL, "FFTW Options", "ex112");
 36:     PetscOptionsEList("-function", "Function type", "ex112", funcNames, NUM_FUNCS, funcNames[function], &func, PETSC_NULL);
 37:     PetscOptionsTruth("-vec_view_draw", "View the functions", "ex112", view, &view, PETSC_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:     /* create vectors of length N=n^DIM */
 49:     for(i = 0, N = 1; i < DIM; i++) N *= dim[i];
 50:     PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n",DIM,N);
 51:     VecCreateSeq(PETSC_COMM_SELF,N,&x);
 52:     PetscObjectSetName((PetscObject) x, "Real space vector");
 53:     VecDuplicate(x,&y);
 54:     PetscObjectSetName((PetscObject) y, "Frequency space vector");
 55:     VecDuplicate(x,&z);
 56:     PetscObjectSetName((PetscObject) z, "Reconstructed vector");
 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;

 64:       VecGetArray(x, &a);
 65:       for(i = 0; i < N; ++i) {
 66:         a[i] = tanh((i - N/2.0)*(10.0/N));
 67:       }
 68:       VecRestoreArray(x, &a);
 69:     }
 70:     if (view) {VecView(x, PETSC_VIEWER_DRAW_WORLD);}

 72:     /* create FFTW object */
 73:     MatCreateSeqFFTW(PETSC_COMM_SELF,DIM,dim,&A);

 75:     /* apply FFTW_FORWARD several times, so the fftw_plan can be reused on different vectors */
 76:     MatMult(A,x,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:       /* apply FFTW_BACKWARD several times */
 81:       MatMultTranspose(A,y,z);
 82:     }
 83: 
 84:     /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
 85:     s = 1.0/(PetscReal)N;
 86:     VecScale(z,s);
 87:     if (view) {VecView(z, PETSC_VIEWER_DRAW_WORLD);}
 88:     VecAXPY(z,-1.0,x);
 89:     VecNorm(z,NORM_1,&enorm);
 90:     if (enorm > 1.e-11){
 91:       PetscPrintf(PETSC_COMM_SELF,"  Error norm of |x - z| %A\n",enorm);
 92:     }

 94:     /* free spaces */
 95:     VecDestroy(x);
 96:     VecDestroy(y);
 97:     VecDestroy(z);
 98:     MatDestroy(A);
 99:   }
100:   PetscRandomDestroy(rdm);
101:   PetscFinalize();
102:   return 0;
103: }