Actual source code: ex121.c
petsc-3.4.2 2013-07-02

      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",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: }