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 {
 16:     RANDOM,
 17:     CONSTANT,
 18:     TANH,
 19:     NUM_FUNCS
 20:   } FuncType;
 21:   const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
 22:   Mat         A, AA;
 23:   PetscMPIInt size;
 24:   PetscInt    N, i, stencil = 1, dof = 1;
 25:   PetscInt    dim[3] = {10, 10, 10}, ndim = 3;
 26:   Vec         coords, x, y, z, xx, yy, zz;
 27:   PetscReal   h[3];
 28:   PetscScalar s;
 29:   PetscRandom rdm;
 30:   PetscReal   norm, enorm;
 31:   PetscInt    func;
 32:   FuncType    function = TANH;
 33:   DM          da, coordsda;
 34:   PetscBool   view_x = PETSC_FALSE, view_y = PETSC_FALSE, view_z = PETSC_FALSE;

 37:   PetscInitialize(&argc, &args, (char *)0, help);
 38:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 40:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "USFFT Options", "ex27");
 41:   PetscOptionsEList("-function", "Function type", "ex27", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
 42:   function = (FuncType)func;
 43:   PetscOptionsEnd();
 44:   PetscOptionsGetBool(NULL, NULL, "-view_x", &view_x, NULL);
 45:   PetscOptionsGetBool(NULL, NULL, "-view_y", &view_y, NULL);
 46:   PetscOptionsGetBool(NULL, NULL, "-view_z", &view_z, NULL);
 47:   PetscOptionsGetIntArray(NULL, NULL, "-dim", dim, &ndim, NULL);

 49:   DMDACreate3d(PETSC_COMM_SELF, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, dim[0], dim[1], dim[2], PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, stencil, NULL, NULL, NULL, &da);
 50:   DMSetFromOptions(da);
 51:   DMSetUp(da);

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

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

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

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

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

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

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

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

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

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

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

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

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

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