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

      1: static char help[] = "Test sequential USFFT interface on a 3-dof field over a uniform DMDA and compares to the result of FFTW acting on a split version of the field\n\n";

      3: /*
      4:   Compiling the code:
      5:       This code uses the complex numbers version of PETSc and the FFTW package, so configure
      6:       must be run to enable these.

      8: */

     10: #define DOF 3

     12: #include <petscmat.h>
     13: #include <petscdmda.h>
     16: PetscInt main(PetscInt argc,char **args)
     17: {
     18:   typedef enum {RANDOM, CONSTANT, TANH, NUM_FUNCS} FuncType;
     19:   const char     *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
     20:   Mat            A, AA;
     21:   PetscMPIInt    size;
     22:   PetscInt       N,i, stencil=1,dof=3;
     23:   PetscInt       dim[3] = {10,10,10}, ndim = 3;
     24:   Vec            coords,x,y,z,xx, yy, zz;
     25:   Vec            xxsplit[DOF], yysplit[DOF], zzsplit[DOF];
     26:   PetscReal      h[3];
     27:   PetscScalar    s;
     28:   PetscRandom    rdm;
     29:   PetscReal      norm, enorm;
     30:   PetscInt       func,ii;
     31:   FuncType       function = TANH;
     32:   DM             da, da1, coordsda;
     33:   PetscBool      view_x = PETSC_FALSE, view_y = PETSC_FALSE, view_z = PETSC_FALSE;

     36:   PetscInitialize(&argc,&args,(char*)0,help);
     37: #if !defined(PETSC_USE_COMPLEX)
     38:   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires complex numbers");
     39: #endif
     40:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
     41:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This is a uniprocessor example only!");
     42:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "USFFT Options", "ex27");
     43:   PetscOptionsEList("-function", "Function type", "ex27", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
     44:   function = (FuncType) func;
     45:   PetscOptionsEnd();
     46:   PetscOptionsGetBool(NULL,"-view_x",&view_x,NULL);
     47:   PetscOptionsGetBool(NULL,"-view_y",&view_y,NULL);
     48:   PetscOptionsGetBool(NULL,"-view_z",&view_z,NULL);
     49:   PetscOptionsGetIntArray(NULL,"-dim",dim,&ndim,NULL);

     51:   /* DMDA with the correct fiber dimension */
     52:   DMDACreate3d(PETSC_COMM_SELF,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,
     53:                       dim[0], dim[1], dim[2],
     54:                       PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
     55:                       dof, stencil,
     56:                       NULL, NULL, NULL,
     57:                       &da);
     58:   /* DMDA with fiber dimension 1 for split fields */
     59:   DMDACreate3d(PETSC_COMM_SELF,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,
     60:                       dim[0], dim[1], dim[2],
     61:                       PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
     62:                       1, stencil,
     63:                       NULL, NULL, NULL,
     64:                       &da1);

     66:   /* Coordinates */
     67:   DMGetCoordinateDM(da, &coordsda);
     68:   DMGetGlobalVector(coordsda, &coords);
     69:   PetscObjectSetName((PetscObject) coords, "Grid coordinates");
     70:   for (i = 0, N = 1; i < 3; i++) {
     71:     h[i] = 1.0/dim[i];
     72:     PetscScalar *a;
     73:     VecGetArray(coords, &a);
     74:     PetscInt j,k,n = 0;
     75:     for (i = 0; i < 3; ++i) {
     76:       for (j = 0; j < dim[i]; ++j) {
     77:         for (k = 0; k < 3; ++k) {
     78:           a[n] = j*h[i]; /* coordinate along the j-th point in the i-th dimension */
     79:           ++n;
     80:         }
     81:       }
     82:     }
     83:     VecRestoreArray(coords, &a);

     85:   }
     86:   DMSetCoordinates(da, coords);
     87:   VecDestroy(&coords);

     89:   /* Work vectors */
     90:   DMGetGlobalVector(da, &x);
     91:   PetscObjectSetName((PetscObject) x, "Real space vector");
     92:   DMGetGlobalVector(da, &xx);
     93:   PetscObjectSetName((PetscObject) xx, "Real space vector");
     94:   DMGetGlobalVector(da, &y);
     95:   PetscObjectSetName((PetscObject) y, "USFFT frequency space vector");
     96:   DMGetGlobalVector(da, &yy);
     97:   PetscObjectSetName((PetscObject) yy, "FFTW frequency space vector");
     98:   DMGetGlobalVector(da, &z);
     99:   PetscObjectSetName((PetscObject) z, "USFFT reconstructed vector");
    100:   DMGetGlobalVector(da, &zz);
    101:   PetscObjectSetName((PetscObject) zz, "FFTW reconstructed vector");
    102:   /* Split vectors for FFTW */
    103:   for (ii = 0; ii < 3; ++ii) {
    104:     DMGetGlobalVector(da1, &xxsplit[ii]);
    105:     PetscObjectSetName((PetscObject) xxsplit[ii], "Real space split vector");
    106:     DMGetGlobalVector(da1, &yysplit[ii]);
    107:     PetscObjectSetName((PetscObject) yysplit[ii], "FFTW frequency space split vector");
    108:     DMGetGlobalVector(da1, &zzsplit[ii]);
    109:     PetscObjectSetName((PetscObject) zzsplit[ii], "FFTW reconstructed split vector");
    110:   }

    113:   PetscPrintf(PETSC_COMM_SELF, "%3-D: USFFT on vector of ");
    114:   for (i = 0, N = 1; i < 3; i++) {
    115:     PetscPrintf(PETSC_COMM_SELF, "dim[%d] = %d ",i,dim[i]);
    116:     N   *= dim[i];
    117:   }
    118:   PetscPrintf(PETSC_COMM_SELF, "; total size %d \n",N);

    121:   if (function == RANDOM) {
    122:     PetscRandomCreate(PETSC_COMM_SELF, &rdm);
    123:     PetscRandomSetFromOptions(rdm);
    124:     VecSetRandom(x, rdm);
    125:     PetscRandomDestroy(&rdm);
    126:   } else if (function == CONSTANT) {
    127:     VecSet(x, 1.0);
    128:   } else if (function == TANH) {
    129:     PetscScalar *a;
    130:     VecGetArray(x, &a);
    131:     PetscInt j,k = 0;
    132:     for (i = 0; i < 3; ++i) {
    133:       for (j = 0; j < dim[i]; ++j) {
    134:         a[k] = tanh((j - dim[i]/2.0)*(10.0/dim[i]));
    135:         ++k;
    136:       }
    137:     }
    138:     VecRestoreArray(x, &a);
    139:   }
    140:   if (view_x) {
    141:     VecView(x, PETSC_VIEWER_STDOUT_WORLD);
    142:   }
    143:   VecCopy(x,xx);
    144:   /* Split xx */
    145:   VecStrideGatherAll(xx,xxsplit, INSERT_VALUES); /*YES! 'Gather' means 'split' (or maybe 'scatter'?)! */

    147:   VecNorm(x,NORM_2,&norm);
    148:   PetscPrintf(PETSC_COMM_SELF, "|x|_2 = %g\n",norm);

    150:   /* create USFFT object */
    151:   MatCreateSeqUSFFT(da,da,&A);
    152:   /* create FFTW object */
    153:   MatCreateSeqFFTW(PETSC_COMM_SELF,3,dim,&AA);

    155:   /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */
    156:   MatMult(A,x,z);
    157:   for (ii = 0; ii < 3; ++ii) {
    158:     MatMult(AA,xxsplit[ii],zzsplit[ii]);
    159:   }
    160:   /* Now apply USFFT and FFTW forward several (3) times */
    161:   for (i=0; i<3; ++i) {
    162:     MatMult(A,x,y);
    163:     for (ii = 0; ii < 3; ++ii) {
    164:       MatMult(AA,xxsplit[ii],yysplit[ii]);
    165:     }
    166:     MatMultTranspose(A,y,z);
    167:     for (ii = 0; ii < 3; ++ii) {
    168:       MatMult(AA,yysplit[ii],zzsplit[ii]);
    169:     }
    170:   }
    171:   /* Unsplit yy */
    172:   VecStrideScatterAll(yysplit, yy, INSERT_VALUES); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */
    173:   /* Unsplit zz */
    174:   VecStrideScatterAll(zzsplit, zz, INSERT_VALUES); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */

    176:   if (view_y) {
    177:     PetscPrintf(PETSC_COMM_WORLD, "y = \n");
    178:     VecView(y, PETSC_VIEWER_STDOUT_WORLD);
    179:     PetscPrintf(PETSC_COMM_WORLD, "yy = \n");
    180:     VecView(yy, PETSC_VIEWER_STDOUT_WORLD);
    181:   }

    183:   if (view_z) {
    184:     PetscPrintf(PETSC_COMM_WORLD, "z = \n");
    185:     VecView(z, PETSC_VIEWER_STDOUT_WORLD);
    186:     PetscPrintf(PETSC_COMM_WORLD, "zz = \n");
    187:     VecView(zz, PETSC_VIEWER_STDOUT_WORLD);
    188:   }

    190:   /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */
    191:   s    = 1.0/(PetscReal)N;
    192:   VecScale(z,s);
    193:   VecAXPY(x,-1.0,z);
    194:   VecNorm(x,NORM_1,&enorm);
    195:   PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n",enorm);

    197:   /* compare xx and zz. FFTW computes an unnormalized DFT, thus zz = N*x */
    198:   s    = 1.0/(PetscReal)N;
    199:   VecScale(zz,s);
    200:   VecAXPY(xx,-1.0,zz);
    201:   VecNorm(xx,NORM_1,&enorm);
    202:   PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n",enorm);

    204:   /* compare y and yy: USFFT and FFTW results*/
    205:   VecNorm(y,NORM_2,&norm);
    206:   VecAXPY(y,-1.0,yy);
    207:   VecNorm(y,NORM_1,&enorm);
    208:   PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n",norm);
    209:   PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n",enorm);

    211:   /* compare z and zz: USFFT and FFTW results*/
    212:   VecNorm(z,NORM_2,&norm);
    213:   VecAXPY(z,-1.0,zz);
    214:   VecNorm(z,NORM_1,&enorm);
    215:   PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n",norm);
    216:   PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n",enorm);

    219:   /* free spaces */
    220:   DMRestoreGlobalVector(da,&x);
    221:   DMRestoreGlobalVector(da,&xx);
    222:   DMRestoreGlobalVector(da,&y);
    223:   DMRestoreGlobalVector(da,&yy);
    224:   DMRestoreGlobalVector(da,&z);
    225:   DMRestoreGlobalVector(da,&zz);

    227:   PetscFinalize();
    228:   return 0;
    229: }