Actual source code: ex121.c

  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;
 24:   FuncType       function = RANDOM;
 25:   PetscTruth     view = PETSC_FALSE;

 28:   PetscInitialize(&argc,&args,(char *)0,help);
 29: #if !defined(PETSC_USE_COMPLEX)
 30:   SETERRQ(PETSC_ERR_SUP, "This example requires complex numbers");
 31: #endif
 32:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 33:   if (size != 1) SETERRQ(PETSC_ERR_SUP, "This is a uniprocessor example only!");
 34:   PetscOptionsBegin(PETSC_COMM_WORLD, PETSC_NULL, "FFTW Options", "ex112");
 35:     PetscOptionsEList("-function", "Function type", "ex112", funcNames, NUM_FUNCS, funcNames[function], &func, PETSC_NULL);
 36:     PetscOptionsTruth("-vec_view_draw", "View the functions", "ex112", view, &view, PETSC_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");
 62:     if (function == RANDOM) {
 63:       VecSetRandom(x, rdm);
 64:     } else if (function == CONSTANT) {
 65:       VecSet(x, 1.0);
 66:     } else if (function == TANH) {
 67:       VecGetArray(x, &a);
 68:       for(i = 0; i < N; ++i) {
 69:         a[i] = tanh((i - N/2.0)*(10.0/N));
 70:       }
 71:       VecRestoreArray(x, &a);
 72:     }
 73:     if (view) {VecView(x, PETSC_VIEWER_DRAW_WORLD);}

 75:     /* Create window function */
 76:     VecGetArray(w, &a);
 77:     for(i = 0; i < N; ++i) {
 78:       // Step Function
 79:       a[i] = (i > N/4 && i < 3*N/4)? 1.0: 0.0;
 80:       // Delta Function
 81:       //a[i] = (i == N/2)? 1.0: 0.0;
 82:     }
 83:     VecRestoreArray(w, &a);
 84:     if (view) {VecView(w, PETSC_VIEWER_DRAW_WORLD);}

 86:     /* create FFTW object */
 87:     //MatCreateSeqFFTW(PETSC_COMM_SELF,DIM,dim,&A);
 88:     MatCreateSeqFFTW(PETSC_COMM_SELF,1,&N,&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);
 96: 
 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 |x - z| %A\n",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: }