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: */

  8:  #include <petscmat.h>

 10: PetscInt main(PetscInt argc,char **args)
 11: {
 12:   typedef enum {RANDOM, CONSTANT, TANH, NUM_FUNCS} FuncType;
 13:   const char     *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
 14:   Mat            A;
 15:   PetscMPIInt    size;
 16:   PetscInt       n = 10,N,ndim=4,dim[4],DIM,i,j;
 17:   Vec            w,x,y1,y2,z1,z2;
 18:   PetscScalar    *a, *a2, *a3;
 19:   PetscScalar    s;
 20:   PetscRandom    rdm;
 21:   PetscReal      enorm;
 22:   PetscInt       func     = 0;
 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", "ex121", 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:     /* create vectors of length N=n^DIM */
 44:     for (i = 0, N = 1; i < DIM; i++) N *= dim[i];
 45:     PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n",DIM,N);
 46:     VecCreateSeq(PETSC_COMM_SELF,N,&x);
 47:     PetscObjectSetName((PetscObject) x, "Real space vector");
 48:     VecDuplicate(x,&w);
 49:     PetscObjectSetName((PetscObject) w, "Window vector");
 50:     VecDuplicate(x,&y1);
 51:     PetscObjectSetName((PetscObject) y1, "Frequency space vector");
 52:     VecDuplicate(x,&y2);
 53:     PetscObjectSetName((PetscObject) y2, "Frequency space window vector");
 54:     VecDuplicate(x,&z1);
 55:     PetscObjectSetName((PetscObject) z1, "Reconstructed convolution");
 56:     VecDuplicate(x,&z2);
 57:     PetscObjectSetName((PetscObject) z2, "Real space convolution");

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

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

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

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

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

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

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

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


140: /*TEST

142:    build:
143:       requires: fftw complex

145:    test:
146:       output_file: output/ex121.out
147:       TODO: Example or FFTW interface is broken

149: TEST*/