Actual source code: ex28.c

  1: static char help[] = "Test sequential USFFT interface on a 3-dof field over a uniform DA 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 petscda.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;
 31:   FuncType       function = TANH;
 32:   DA             da, da1, coordsda;
 33:   PetscTruth     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_ERR_SUP, "This example requires complex numbers");
 39: #endif
 40:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 41:   if (size != 1) SETERRQ(PETSC_ERR_SUP, "This is a uniprocessor example only!");
 42:   PetscOptionsBegin(PETSC_COMM_WORLD, PETSC_NULL, "USFFT Options", "ex27");
 43:     PetscOptionsEList("-function", "Function type", "ex27", funcNames, NUM_FUNCS, funcNames[function], &func, PETSC_NULL);
 44:     function = (FuncType) func;
 45:   PetscOptionsEnd();
 46:   PetscOptionsGetTruth(PETSC_NULL,"-view_x",&view_x,PETSC_NULL);
 47:   PetscOptionsGetTruth(PETSC_NULL,"-view_y",&view_y,PETSC_NULL);
 48:   PetscOptionsGetTruth(PETSC_NULL,"-view_z",&view_z,PETSC_NULL);
 49:   PetscOptionsGetIntArray(PETSC_NULL,"-dim",dim,&ndim,PETSC_NULL);

 51:   // DA with the correct fiber dimension
 52:   DACreate3d(PETSC_COMM_SELF,DA_NONPERIODIC,DA_STENCIL_STAR,
 53:                     dim[0], dim[1], dim[2],
 54:                     PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
 55:                     dof, stencil,
 56:                     PETSC_NULL, PETSC_NULL, PETSC_NULL,
 57:                     &da);
 58: 
 59:   // DA with fiber dimension 1 for split fields
 60:   DACreate3d(PETSC_COMM_SELF,DA_NONPERIODIC,DA_STENCIL_STAR,
 61:                     dim[0], dim[1], dim[2],
 62:                     PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
 63:                     1, stencil,
 64:                     PETSC_NULL, PETSC_NULL, PETSC_NULL,
 65:                     &da1);
 66: 
 67: 
 68:   // Coordinates
 69:   DAGetCoordinateDA(da, &coordsda);
 70:   DAGetGlobalVector(coordsda, &coords);
 71:   PetscObjectSetName((PetscObject) coords, "Grid coordinates");
 72:   for(i = 0, N = 1; i < 3; i++) {
 73:     h[i] = 1.0/dim[i];
 74:     PetscScalar *a;
 75:     VecGetArray(coords, &a);
 76:     PetscInt j,k,n = 0;
 77:     for(i = 0; i < 3; ++i) {
 78:       for(j = 0; j < dim[i]; ++j){
 79:         for(k = 0; k < 3; ++k) {
 80:           a[n] = j*h[i]; // coordinate along the j-th point in the i-th dimension
 81:           ++n;
 82:         }
 83:       }
 84:     }
 85:     VecRestoreArray(coords, &a);

 87:   }
 88:   DASetCoordinates(da, coords);
 89:   VecDestroy(coords);

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


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

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

151:   VecNorm(x,NORM_2,&norm);
152:   PetscPrintf(PETSC_COMM_SELF, "|x|_2 = %g\n",norm);
153: 
154:   /* create USFFT object */
155:   MatCreateSeqUSFFT(da,da,&A);
156:   /* create FFTW object */
157:   MatCreateSeqFFTW(PETSC_COMM_SELF,3,dim,&AA);
158: 
159:   /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */
160:   MatMult(A,x,z);
161:   for(int ii = 0; ii < 3; ++ii) {
162:     MatMult(AA,xxsplit[ii],zzsplit[ii]);
163:   }
164:   // Now apply USFFT and FFTW forward several (3) times
165:   for (i=0; i<3; ++i){
166:     MatMult(A,x,y);
167:     for(int ii = 0; ii < 3; ++ii) {
168:       MatMult(AA,xxsplit[ii],yysplit[ii]);
169:     }
170:     MatMultTranspose(A,y,z);
171:     for(int ii = 0; ii < 3; ++ii) {
172:       MatMult(AA,yysplit[ii],zzsplit[ii]);
173:     }
174:   }
175:   // Unsplit yy
176:   VecStrideScatterAll(yysplit, yy, INSERT_VALUES); //YES! 'Scatter' means 'collect' (or maybe 'gather'?)!
177:   // Unsplit zz
178:   VecStrideScatterAll(zzsplit, zz, INSERT_VALUES); //YES! 'Scatter' means 'collect' (or maybe 'gather'?)!

180:   if(view_y) {
181:     PetscPrintf(PETSC_COMM_WORLD, "y = \n");
182:     VecView(y, PETSC_VIEWER_STDOUT_WORLD);
183:     PetscPrintf(PETSC_COMM_WORLD, "yy = \n");
184:     VecView(yy, PETSC_VIEWER_STDOUT_WORLD);
185:   }
186: 
187:   if(view_z) {
188:     PetscPrintf(PETSC_COMM_WORLD, "z = \n");
189:     VecView(z, PETSC_VIEWER_STDOUT_WORLD);
190:     PetscPrintf(PETSC_COMM_WORLD, "zz = \n");
191:     VecView(zz, PETSC_VIEWER_STDOUT_WORLD);
192:   }
193: 
194:   /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */
195:   s = 1.0/(PetscReal)N;
196:   VecScale(z,s);
197:   VecAXPY(x,-1.0,z);
198:   VecNorm(x,NORM_1,&enorm);
199:   PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n",enorm);

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

208:   /* compare y and yy: USFFT and FFTW results*/
209:   VecNorm(y,NORM_2,&norm);
210:   VecAXPY(y,-1.0,yy);
211:   VecNorm(y,NORM_1,&enorm);
212:   PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n",norm);
213:   PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n",enorm);
214: 
215:   /* compare z and zz: USFFT and FFTW results*/
216:   VecNorm(z,NORM_2,&norm);
217:   VecAXPY(z,-1.0,zz);
218:   VecNorm(z,NORM_1,&enorm);
219:   PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n",norm);
220:   PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n",enorm);
221: 

223:   /* free spaces */
224:   DARestoreGlobalVector(da,&x);
225:   DARestoreGlobalVector(da,&xx);
226:   DARestoreGlobalVector(da,&y);
227:   DARestoreGlobalVector(da,&yy);
228:   DARestoreGlobalVector(da,&z);
229:   DARestoreGlobalVector(da,&zz);

231:   PetscFinalize();
232:   return 0;
233: }