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>

 11: PetscInt main(PetscInt 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,j;
 18:   Vec            w,x,y1,y2,z1,z2;
 19:   PetscScalar    *a, *a2, *a3;
 20:   PetscScalar    s;
 21:   PetscRandom    rdm;
 22:   PetscReal      enorm;
 23:   PetscInt       func     = 0;
 24:   FuncType       function = RANDOM;
 25:   PetscBool      view     = PETSC_FALSE;

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

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

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

 63:     if (function == RANDOM) {
 64:       VecSetRandom(x, rdm);
 65:     } else if (function == CONSTANT) {
 66:       VecSet(x, 1.0);
 67:     } else if (function == TANH) {
 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:     /* Create window function */
 77:     VecGetArray(w, &a);
 78:     for (i = 0; i < N; ++i) {
 79:       /* Step Function */
 80:       a[i] = (i > N/4 && i < 3*N/4) ? 1.0 : 0.0;
 81:       /* Delta Function */
 82:       /*a[i] = (i == N/2)? 1.0: 0.0; */
 83:     }
 84:     VecRestoreArray(w, &a);
 85:     if (view) {VecView(w, PETSC_VIEWER_DRAW_WORLD);}

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

 90:     /* Convolve x with w*/
 91:     MatMult(A,x,y1);
 92:     MatMult(A,w,y2);
 93:     VecPointwiseMult(y1, y1, y2);
 94:     if (view && i == 0) {VecView(y1, PETSC_VIEWER_DRAW_WORLD);}
 95:     MatMultTranspose(A,y1,z1);

 97:     /* Compute the real space convolution */
 98:     VecGetArray(x, &a);
 99:     VecGetArray(w, &a2);
100:     VecGetArray(z2, &a3);
101:     for (i = 0; i < N; ++i) {
102:       /* PetscInt checkInd = (i > N/2-1)? i-N/2: i+N/2;*/

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

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

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

129:     /* free spaces */
130:     VecDestroy(&x);
131:     VecDestroy(&y1);
132:     VecDestroy(&y2);
133:     VecDestroy(&z1);
134:     VecDestroy(&z2);
135:     VecDestroy(&w);
136:     MatDestroy(&A);
137:   }
138:   PetscRandomDestroy(&rdm);
139:   PetscFinalize();
140:   return 0;
141: }