Actual source code: ex28.c

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

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

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

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

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

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


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


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

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

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

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

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

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

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

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

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

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


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

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