1: static char help[] = "Test sequential FFTW convolution\n\n";

  3: /*
  4:   Compiling the code:
  5:     This code uses the complex numbers, so configure must be given --with-scalar-type=complex to enable this
  6: */
  7:  #include <petscmat.h>

  9: PetscInt main(PetscInt argc,char **args)
 10: {
 11:   typedef enum {RANDOM, CONSTANT, TANH, NUM_FUNCS} FuncType;
 12:   const char     *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
 13:   Mat            A;
 14:   PetscMPIInt    size;
 15:   PetscInt       n = 10,N,ndim=4,dim[4],DIM,i,j;
 16:   Vec            w,x,y1,y2,z1,z2;
 17:   PetscScalar    *a, *a2, *a3;
 18:   PetscScalar    s;
 19:   PetscRandom    rdm;
 20:   PetscReal      enorm;
 21:   PetscInt       func     = 0;
 22:   FuncType       function = RANDOM;
 23:   PetscBool      view     = PETSC_FALSE;

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

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

 41:   for (DIM = 1; DIM < 5; DIM++) {
 42:     /* create vectors of length N=n^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);
 45:     VecCreateSeq(PETSC_COMM_SELF,N,&x);
 46:     PetscObjectSetName((PetscObject) x, "Real space vector");
 47:     VecDuplicate(x,&w);
 48:     PetscObjectSetName((PetscObject) w, "Window vector");
 49:     VecDuplicate(x,&y1);
 50:     PetscObjectSetName((PetscObject) y1, "Frequency space vector");
 51:     VecDuplicate(x,&y2);
 52:     PetscObjectSetName((PetscObject) y2, "Frequency space window vector");
 53:     VecDuplicate(x,&z1);
 54:     PetscObjectSetName((PetscObject) z1, "Reconstructed convolution");
 55:     VecDuplicate(x,&z2);
 56:     PetscObjectSetName((PetscObject) z2, "Real space convolution");

 58:     if (function == RANDOM) {
 59:       VecSetRandom(x, rdm);
 60:     } else if (function == CONSTANT) {
 61:       VecSet(x, 1.0);
 62:     } else if (function == TANH) {
 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:     /* Create window function */
 72:     VecGetArray(w, &a);
 73:     for (i = 0; i < N; ++i) {
 74:       /* Step Function */
 75:       a[i] = (i > N/4 && i < 3*N/4) ? 1.0 : 0.0;
 76:       /* Delta Function */
 77:       /*a[i] = (i == N/2)? 1.0: 0.0; */
 78:     }
 79:     VecRestoreArray(w, &a);
 80:     if (view) {VecView(w, PETSC_VIEWER_DRAW_WORLD);}

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

 85:     /* Convolve x with w*/
 86:     MatMult(A,x,y1);
 87:     MatMult(A,w,y2);
 88:     VecPointwiseMult(y1, y1, y2);
 89:     if (view && i == 0) {VecView(y1, PETSC_VIEWER_DRAW_WORLD);}
 90:     MatMultTranspose(A,y1,z1);

 92:     /* Compute the real space convolution */
 93:     VecGetArray(x, &a);
 94:     VecGetArray(w, &a2);
 95:     VecGetArray(z2, &a3);
 96:     for (i = 0; i < N; ++i) {
 97:       /* PetscInt checkInd = (i > N/2-1)? i-N/2: i+N/2;*/

 99:       /*if (!(i%100)) PetscPrintf(PETSC_COMM_WORLD, "Finished convolution row %d\n", i);*/
100:       a3[i] = 0.0;
101:       for (j = -N/2+1; j < N/2; ++j) {
102:         PetscInt xpInd   = (j < 0) ? N+j : j;
103:         PetscInt diffInd = (i-j < 0) ? N-(j-i) : (i-j > N-1) ? i-j-N : i-j;

105:         a3[i] += a[xpInd]*a2[diffInd];
106:       }
107:       /*if (PetscAbsScalar(a3[i]) > PetscAbsScalar(a[checkInd])+0.1) PetscPrintf(PETSC_COMM_WORLD, "Invalid convolution at row %d\n", i);*/
108:     }
109:     VecRestoreArray(x, &a);
110:     VecRestoreArray(w, &a2);
111:     VecRestoreArray(z2, &a3);

113:     /* compare z1 and z2. FFTW computes an unnormalized DFT, thus z1 = N*z2 */
114:     s    = 1.0/(PetscReal)N;
115:     VecScale(z1,s);
116:     if (view) {VecView(z1, PETSC_VIEWER_DRAW_WORLD);}
117:     if (view) {VecView(z2, PETSC_VIEWER_DRAW_WORLD);}
118:     VecAXPY(z1,-1.0,z2);
119:     VecNorm(z1,NORM_1,&enorm);
120:     if (enorm > 1.e-11) {
121:       PetscPrintf(PETSC_COMM_SELF,"  Error norm of |z1 - z2| %g\n",(double)enorm);
122:     }

124:     /* free spaces */
125:     VecDestroy(&x);
126:     VecDestroy(&y1);
127:     VecDestroy(&y2);
128:     VecDestroy(&z1);
129:     VecDestroy(&z2);
130:     VecDestroy(&w);
131:     MatDestroy(&A);
132:   }
133:   PetscRandomDestroy(&rdm);
134:   PetscFinalize();
135:   return ierr;
136: }