Actual source code: ex27.c

  1: static char help[] = "Test sequential USFFT interface on a uniform DMDA and compares the result to FFTW\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: #include <petscmat.h>
 11: #include <petscdm.h>
 12: #include <petscdmda.h>
 13: int main(int argc,char **args)
 14: {
 15:   typedef enum {RANDOM, CONSTANT, TANH, NUM_FUNCS} FuncType;
 16:   const char     *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
 17:   Mat            A, AA;
 18:   PetscMPIInt    size;
 19:   PetscInt       N,i, stencil=1,dof=1;
 20:   PetscInt       dim[3] = {10,10,10}, ndim = 3;
 21:   Vec            coords,x,y,z,xx,yy,zz;
 22:   PetscReal      h[3];
 23:   PetscScalar    s;
 24:   PetscRandom    rdm;
 25:   PetscReal      norm, enorm;
 26:   PetscInt       func;
 27:   FuncType       function = TANH;
 28:   DM             da, coordsda;
 29:   PetscBool      view_x = PETSC_FALSE, view_y = PETSC_FALSE, view_z = PETSC_FALSE;

 32:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 33:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 34:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This is a uniprocessor example only!");
 35:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "USFFT Options", "ex27");
 36:   PetscOptionsEList("-function", "Function type", "ex27", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
 37:   function = (FuncType) func;
 38:   PetscOptionsEnd();
 39:   PetscOptionsGetBool(NULL,NULL,"-view_x",&view_x,NULL);
 40:   PetscOptionsGetBool(NULL,NULL,"-view_y",&view_y,NULL);
 41:   PetscOptionsGetBool(NULL,NULL,"-view_z",&view_z,NULL);
 42:   PetscOptionsGetIntArray(NULL,NULL,"-dim",dim,&ndim,NULL);

 44:   DMDACreate3d(PETSC_COMM_SELF,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,dim[0], dim[1], dim[2],
 45:                       PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,dof, stencil,NULL, NULL, NULL,&da);
 46:   DMSetFromOptions(da);
 47:   DMSetUp(da);

 49:   /* Coordinates */
 50:   DMGetCoordinateDM(da, &coordsda);
 51:   DMGetGlobalVector(coordsda, &coords);
 52:   PetscObjectSetName((PetscObject) coords, "Grid coordinates");
 53:   for (i = 0, N = 1; i < 3; i++) {
 54:     h[i] = 1.0/dim[i];
 55:     PetscScalar *a;
 56:     VecGetArray(coords, &a);
 57:     PetscInt j,k,n = 0;
 58:     for (i = 0; i < 3; ++i) {
 59:       for (j = 0; j < dim[i]; ++j) {
 60:         for (k = 0; k < 3; ++k) {
 61:           a[n] = j*h[i]; /* coordinate along the j-th point in the i-th dimension */
 62:           ++n;
 63:         }
 64:       }
 65:     }
 66:     VecRestoreArray(coords, &a);

 68:   }
 69:   DMSetCoordinates(da, coords);

 71:   /* Work vectors */
 72:   DMGetGlobalVector(da, &x);
 73:   PetscObjectSetName((PetscObject) x, "Real space vector");
 74:   DMGetGlobalVector(da, &xx);
 75:   PetscObjectSetName((PetscObject) xx, "Real space vector");
 76:   DMGetGlobalVector(da, &y);
 77:   PetscObjectSetName((PetscObject) y, "USFFT frequency space vector");
 78:   DMGetGlobalVector(da, &yy);
 79:   PetscObjectSetName((PetscObject) yy, "FFTW frequency space vector");
 80:   DMGetGlobalVector(da, &z);
 81:   PetscObjectSetName((PetscObject) z, "USFFT reconstructed vector");
 82:   DMGetGlobalVector(da, &zz);
 83:   PetscObjectSetName((PetscObject) zz, "FFTW reconstructed vector");

 85:   PetscPrintf(PETSC_COMM_SELF, "%3-D: USFFT on vector of ");
 86:   for (i = 0, N = 1; i < 3; i++) {
 87:     PetscPrintf(PETSC_COMM_SELF, "dim[%d] = %d ",i,dim[i]);
 88:     N   *= dim[i];
 89:   }
 90:   PetscPrintf(PETSC_COMM_SELF, "; total size %d \n",N);

 92:   if (function == RANDOM) {
 93:     PetscRandomCreate(PETSC_COMM_SELF, &rdm);
 94:     PetscRandomSetFromOptions(rdm);
 95:     VecSetRandom(x, rdm);
 96:     PetscRandomDestroy(&rdm);
 97:   } else if (function == CONSTANT) {
 98:     VecSet(x, 1.0);
 99:   } else if (function == TANH) {
100:     PetscScalar *a;
101:     VecGetArray(x, &a);
102:     PetscInt j,k = 0;
103:     for (i = 0; i < 3; ++i) {
104:       for (j = 0; j < dim[i]; ++j) {
105:         a[k] = tanh((j - dim[i]/2.0)*(10.0/dim[i]));
106:         ++k;
107:       }
108:     }
109:     VecRestoreArray(x, &a);
110:   }
111:   if (view_x) {
112:     VecView(x, PETSC_VIEWER_STDOUT_WORLD);
113:   }
114:   VecCopy(x,xx);

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

119:   /* create USFFT object */
120:   MatCreateSeqUSFFT(coords,da,&A);
121:   /* create FFTW object */
122:   MatCreateSeqFFTW(PETSC_COMM_SELF,3,dim,&AA);

124:   /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */
125:   MatMult(A,x,z);
126:   MatMult(AA,xx,zz);
127:   /* Now apply USFFT and FFTW forward several (3) times */
128:   for (i=0; i<3; ++i) {
129:     MatMult(A,x,y);
130:     MatMult(AA,xx,yy);
131:     MatMultTranspose(A,y,z);
132:     MatMultTranspose(AA,yy,zz);
133:   }

135:   if (view_y) {
136:     PetscPrintf(PETSC_COMM_WORLD, "y = \n");
137:     VecView(y, PETSC_VIEWER_STDOUT_WORLD);
138:     PetscPrintf(PETSC_COMM_WORLD, "yy = \n");
139:     VecView(yy, PETSC_VIEWER_STDOUT_WORLD);
140:   }

142:   if (view_z) {
143:     PetscPrintf(PETSC_COMM_WORLD, "z = \n");
144:     VecView(z, PETSC_VIEWER_STDOUT_WORLD);
145:     PetscPrintf(PETSC_COMM_WORLD, "zz = \n");
146:     VecView(zz, PETSC_VIEWER_STDOUT_WORLD);
147:   }

149:   /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */
150:   s    = 1.0/(PetscReal)N;
151:   VecScale(z,s);
152:   VecAXPY(x,-1.0,z);
153:   VecNorm(x,NORM_1,&enorm);
154:   PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n",enorm);

156:   /* compare xx and zz. FFTW computes an unnormalized DFT, thus zz = N*x */
157:   s    = 1.0/(PetscReal)N;
158:   VecScale(zz,s);
159:   VecAXPY(xx,-1.0,zz);
160:   VecNorm(xx,NORM_1,&enorm);
161:   PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n",enorm);

163:   /* compare y and yy: USFFT and FFTW results*/
164:   VecNorm(y,NORM_2,&norm);
165:   VecAXPY(y,-1.0,yy);
166:   VecNorm(y,NORM_1,&enorm);
167:   PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n",norm);
168:   PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n",enorm);

170:   /* compare z and zz: USFFT and FFTW results*/
171:   VecNorm(z,NORM_2,&norm);
172:   VecAXPY(z,-1.0,zz);
173:   VecNorm(z,NORM_1,&enorm);
174:   PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n",norm);
175:   PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n",enorm);

177:   /* free spaces */
178:   DMRestoreGlobalVector(da,&x);
179:   DMRestoreGlobalVector(da,&xx);
180:   DMRestoreGlobalVector(da,&y);
181:   DMRestoreGlobalVector(da,&yy);
182:   DMRestoreGlobalVector(da,&z);
183:   DMRestoreGlobalVector(da,&zz);
184:   VecDestroy(&coords);
185:   DMDestroy(&da);
186:   PetscFinalize();
187:   return ierr;
188: }