Actual source code: ex27.c

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

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

 48: 

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

 75:   }
 76:   DASetCoordinates(da, coords);

 78:   // Work vectors
 79:   DAGetGlobalVector(da, &x);
 80:   PetscObjectSetName((PetscObject) x, "Real space vector");
 81:   DAGetGlobalVector(da, &xx);
 82:   PetscObjectSetName((PetscObject) xx, "Real space vector");
 83:   DAGetGlobalVector(da, &y);
 84:   PetscObjectSetName((PetscObject) y, "USFFT frequency space vector");
 85:   DAGetGlobalVector(da, &yy);
 86:   PetscObjectSetName((PetscObject) yy, "FFTW frequency space vector");
 87:   DAGetGlobalVector(da, &z);
 88:   PetscObjectSetName((PetscObject) z, "USFFT reconstructed vector");
 89:   DAGetGlobalVector(da, &zz);
 90:   PetscObjectSetName((PetscObject) zz, "FFTW reconstructed vector");

 92:   PetscPrintf(PETSC_COMM_SELF, "%3-D: USFFT on vector of ");
 93:   for(i = 0, N = 1; i < 3; i++) {
 94:     PetscPrintf(PETSC_COMM_SELF, "dim[%d] = %d ",i,dim[i]);
 95:     N *= dim[i];
 96:   }
 97:   PetscPrintf(PETSC_COMM_SELF, "; total size %d \n",N);

 99: 
100:   if (function == RANDOM) {
101:     PetscRandomCreate(PETSC_COMM_SELF, &rdm);
102:     PetscRandomSetFromOptions(rdm);
103:     VecSetRandom(x, rdm);
104:     PetscRandomDestroy(rdm);
105:   }
106:   else if (function == CONSTANT) {
107:     VecSet(x, 1.0);
108:   }
109:   else if (function == TANH) {
110:     PetscScalar *a;
111:     VecGetArray(x, &a);
112:     PetscInt j,k = 0;
113:     for(i = 0; i < 3; ++i) {
114:       for(j = 0; j < dim[i]; ++j) {
115:         a[k] = tanh((j - dim[i]/2.0)*(10.0/dim[i]));
116:         ++k;
117:       }
118:     }
119:     VecRestoreArray(x, &a);
120:   }
121:   if(view_x) {
122:     VecView(x, PETSC_VIEWER_STDOUT_WORLD);
123:   }
124:   VecCopy(x,xx);

126:   VecNorm(x,NORM_2,&norm);
127:   PetscPrintf(PETSC_COMM_SELF, "|x|_2 = %g\n",norm);
128: 
129:   /* create USFFT object */
130:   MatCreateSeqUSFFT(coords,da,&A);
131:   /* create FFTW object */
132:   MatCreateSeqFFTW(PETSC_COMM_SELF,3,dim,&AA);
133: 
134:   /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */
135:   MatMult(A,x,z);
136:   MatMult(AA,xx,zz);
137:   // Now apply USFFT and FFTW forward several (3) times
138:   for (i=0; i<3; ++i){
139:     MatMult(A,x,y);
140:     MatMult(AA,xx,yy);
141:     MatMultTranspose(A,y,z);
142:     MatMultTranspose(AA,yy,zz);
143:   }

145:   if(view_y) {
146:     PetscPrintf(PETSC_COMM_WORLD, "y = \n");
147:     VecView(y, PETSC_VIEWER_STDOUT_WORLD);
148:     PetscPrintf(PETSC_COMM_WORLD, "yy = \n");
149:     VecView(yy, PETSC_VIEWER_STDOUT_WORLD);
150:   }
151: 
152:   if(view_z) {
153:     PetscPrintf(PETSC_COMM_WORLD, "z = \n");
154:     VecView(z, PETSC_VIEWER_STDOUT_WORLD);
155:     PetscPrintf(PETSC_COMM_WORLD, "zz = \n");
156:     VecView(zz, PETSC_VIEWER_STDOUT_WORLD);
157:   }
158: 
159:   /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */
160:   s = 1.0/(PetscReal)N;
161:   VecScale(z,s);
162:   VecAXPY(x,-1.0,z);
163:   VecNorm(x,NORM_1,&enorm);
164:   PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n",enorm);

166:   /* compare xx and zz. FFTW computes an unnormalized DFT, thus zz = N*x */
167:   s = 1.0/(PetscReal)N;
168:   VecScale(zz,s);
169:   VecAXPY(xx,-1.0,zz);
170:   VecNorm(xx,NORM_1,&enorm);
171:   PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n",enorm);

173:   /* compare y and yy: USFFT and FFTW results*/
174:   VecNorm(y,NORM_2,&norm);
175:   VecAXPY(y,-1.0,yy);
176:   VecNorm(y,NORM_1,&enorm);
177:   PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n",norm);
178:   PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n",enorm);
179: 
180:   /* compare z and zz: USFFT and FFTW results*/
181:   VecNorm(z,NORM_2,&norm);
182:   VecAXPY(z,-1.0,zz);
183:   VecNorm(z,NORM_1,&enorm);
184:   PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n",norm);
185:   PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n",enorm);
186: 

188:   /* free spaces */
189:   DARestoreGlobalVector(da,&x);
190:   DARestoreGlobalVector(da,&xx);
191:   DARestoreGlobalVector(da,&y);
192:   DARestoreGlobalVector(da,&yy);
193:   DARestoreGlobalVector(da,&z);
194:   DARestoreGlobalVector(da,&zz);
195:   VecDestroy(coords);
196:   DADestroy(da);
197:   PetscFinalize();
198:   return 0;
199: }