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;

 36:   PetscFunctionBeginUser;
 37:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 38:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 39:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only!");
 40:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "USFFT Options", "ex27");
 41:   PetscCall(PetscOptionsEList("-function", "Function type", "ex27", funcNames, NUM_FUNCS, funcNames[function], &func, NULL));
 42:   function = (FuncType)func;
 43:   PetscOptionsEnd();
 44:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_x", &view_x, NULL));
 45:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_y", &view_y, NULL));
 46:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_z", &view_z, NULL));
 47:   PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dim", dim, &ndim, NULL));

 49:   PetscCall(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:   PetscCall(DMSetFromOptions(da));
 51:   PetscCall(DMSetUp(da));

 53:   /* Coordinates */
 54:   PetscCall(DMGetCoordinateDM(da, &coordsda));
 55:   PetscCall(DMGetGlobalVector(coordsda, &coords));
 56:   PetscCall(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:     PetscCall(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:     PetscCall(VecRestoreArray(coords, &a));
 71:   }
 72:   PetscCall(DMSetCoordinates(da, coords));

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

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

 95:   if (function == RANDOM) {
 96:     PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
 97:     PetscCall(PetscRandomSetFromOptions(rdm));
 98:     PetscCall(VecSetRandom(x, rdm));
 99:     PetscCall(PetscRandomDestroy(&rdm));
100:   } else if (function == CONSTANT) {
101:     PetscCall(VecSet(x, 1.0));
102:   } else if (function == TANH) {
103:     PetscScalar *a;
104:     PetscCall(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:     PetscCall(VecRestoreArray(x, &a));
113:   }
114:   if (view_x) PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
115:   PetscCall(VecCopy(x, xx));

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

120:   /* create USFFT object */
121:   PetscCall(MatCreateSeqUSFFT(coords, da, &A));
122:   /* create FFTW object */
123:   PetscCall(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:   PetscCall(MatMult(A, x, z));
127:   PetscCall(MatMult(AA, xx, zz));
128:   /* Now apply USFFT and FFTW forward several (3) times */
129:   for (i = 0; i < 3; ++i) {
130:     PetscCall(MatMult(A, x, y));
131:     PetscCall(MatMult(AA, xx, yy));
132:     PetscCall(MatMultTranspose(A, y, z));
133:     PetscCall(MatMultTranspose(AA, yy, zz));
134:   }

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

143:   if (view_z) {
144:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "z = \n"));
145:     PetscCall(VecView(z, PETSC_VIEWER_STDOUT_WORLD));
146:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "zz = \n"));
147:     PetscCall(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:   PetscCall(VecScale(z, s));
153:   PetscCall(VecAXPY(x, -1.0, z));
154:   PetscCall(VecNorm(x, NORM_1, &enorm));
155:   PetscCall(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:   PetscCall(VecScale(zz, s));
160:   PetscCall(VecAXPY(xx, -1.0, zz));
161:   PetscCall(VecNorm(xx, NORM_1, &enorm));
162:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n", enorm));

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

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

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