Actual source code: ex28.c

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

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

 52:   /* DMDA with the correct fiber dimension */
 53:   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));
 54:   PetscCall(DMSetFromOptions(da));
 55:   PetscCall(DMSetUp(da));
 56:   /* DMDA with fiber dimension 1 for split fields */
 57:   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, 1, stencil, NULL, NULL, NULL, &da1));
 58:   PetscCall(DMSetFromOptions(da1));
 59:   PetscCall(DMSetUp(da1));

 61:   /* Coordinates */
 62:   PetscCall(DMGetCoordinateDM(da, &coordsda));
 63:   PetscCall(DMGetGlobalVector(coordsda, &coords));
 64:   PetscCall(PetscObjectSetName((PetscObject)coords, "Grid coordinates"));
 65:   for (i = 0, N = 1; i < 3; i++) {
 66:     h[i] = 1.0 / dim[i];
 67:     PetscScalar *a;
 68:     PetscCall(VecGetArray(coords, &a));
 69:     PetscInt j, k, n = 0;
 70:     for (i = 0; i < 3; ++i) {
 71:       for (j = 0; j < dim[i]; ++j) {
 72:         for (k = 0; k < 3; ++k) {
 73:           a[n] = j * h[i]; /* coordinate along the j-th point in the i-th dimension */
 74:           ++n;
 75:         }
 76:       }
 77:     }
 78:     PetscCall(VecRestoreArray(coords, &a));
 79:   }
 80:   PetscCall(DMSetCoordinates(da, coords));
 81:   PetscCall(VecDestroy(&coords));

 83:   /* Work vectors */
 84:   PetscCall(DMGetGlobalVector(da, &x));
 85:   PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector"));
 86:   PetscCall(DMGetGlobalVector(da, &xx));
 87:   PetscCall(PetscObjectSetName((PetscObject)xx, "Real space vector"));
 88:   PetscCall(DMGetGlobalVector(da, &y));
 89:   PetscCall(PetscObjectSetName((PetscObject)y, "USFFT frequency space vector"));
 90:   PetscCall(DMGetGlobalVector(da, &yy));
 91:   PetscCall(PetscObjectSetName((PetscObject)yy, "FFTW frequency space vector"));
 92:   PetscCall(DMGetGlobalVector(da, &z));
 93:   PetscCall(PetscObjectSetName((PetscObject)z, "USFFT reconstructed vector"));
 94:   PetscCall(DMGetGlobalVector(da, &zz));
 95:   PetscCall(PetscObjectSetName((PetscObject)zz, "FFTW reconstructed vector"));
 96:   /* Split vectors for FFTW */
 97:   for (ii = 0; ii < 3; ++ii) {
 98:     PetscCall(DMGetGlobalVector(da1, &xxsplit[ii]));
 99:     PetscCall(PetscObjectSetName((PetscObject)xxsplit[ii], "Real space split vector"));
100:     PetscCall(DMGetGlobalVector(da1, &yysplit[ii]));
101:     PetscCall(PetscObjectSetName((PetscObject)yysplit[ii], "FFTW frequency space split vector"));
102:     PetscCall(DMGetGlobalVector(da1, &zzsplit[ii]));
103:     PetscCall(PetscObjectSetName((PetscObject)zzsplit[ii], "FFTW reconstructed split vector"));
104:   }

106:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "%3-" PetscInt_FMT ": USFFT on vector of "));
107:   for (i = 0, N = 1; i < 3; i++) {
108:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "dim[%d] = %d ", i, dim[i]));
109:     N *= dim[i];
110:   }
111:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "; total size %d \n", N));

113:   if (function == RANDOM) {
114:     PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
115:     PetscCall(PetscRandomSetFromOptions(rdm));
116:     PetscCall(VecSetRandom(x, rdm));
117:     PetscCall(PetscRandomDestroy(&rdm));
118:   } else if (function == CONSTANT) {
119:     PetscCall(VecSet(x, 1.0));
120:   } else if (function == TANH) {
121:     PetscScalar *a;
122:     PetscCall(VecGetArray(x, &a));
123:     PetscInt j, k = 0;
124:     for (i = 0; i < 3; ++i) {
125:       for (j = 0; j < dim[i]; ++j) {
126:         a[k] = tanh((j - dim[i] / 2.0) * (10.0 / dim[i]));
127:         ++k;
128:       }
129:     }
130:     PetscCall(VecRestoreArray(x, &a));
131:   }
132:   if (view_x) PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
133:   PetscCall(VecCopy(x, xx));
134:   /* Split xx */
135:   PetscCall(VecStrideGatherAll(xx, xxsplit, INSERT_VALUES)); /*YES! 'Gather' means 'split' (or maybe 'scatter'?)! */

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

140:   /* create USFFT object */
141:   PetscCall(MatCreateSeqUSFFT(da, da, &A));
142:   /* create FFTW object */
143:   PetscCall(MatCreateSeqFFTW(PETSC_COMM_SELF, 3, dim, &AA));

145:   /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */
146:   PetscCall(MatMult(A, x, z));
147:   for (ii = 0; ii < 3; ++ii) PetscCall(MatMult(AA, xxsplit[ii], zzsplit[ii]));
148:   /* Now apply USFFT and FFTW forward several (3) times */
149:   for (i = 0; i < 3; ++i) {
150:     PetscCall(MatMult(A, x, y));
151:     for (ii = 0; ii < 3; ++ii) PetscCall(MatMult(AA, xxsplit[ii], yysplit[ii]));
152:     PetscCall(MatMultTranspose(A, y, z));
153:     for (ii = 0; ii < 3; ++ii) PetscCall(MatMult(AA, yysplit[ii], zzsplit[ii]));
154:   }
155:   /* Unsplit yy */
156:   PetscCall(VecStrideScatterAll(yysplit, yy, INSERT_VALUES)); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */
157:   /* Unsplit zz */
158:   PetscCall(VecStrideScatterAll(zzsplit, zz, INSERT_VALUES)); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */

160:   if (view_y) {
161:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "y = \n"));
162:     PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD));
163:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "yy = \n"));
164:     PetscCall(VecView(yy, PETSC_VIEWER_STDOUT_WORLD));
165:   }

167:   if (view_z) {
168:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "z = \n"));
169:     PetscCall(VecView(z, PETSC_VIEWER_STDOUT_WORLD));
170:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "zz = \n"));
171:     PetscCall(VecView(zz, PETSC_VIEWER_STDOUT_WORLD));
172:   }

174:   /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */
175:   s = 1.0 / (PetscReal)N;
176:   PetscCall(VecScale(z, s));
177:   PetscCall(VecAXPY(x, -1.0, z));
178:   PetscCall(VecNorm(x, NORM_1, &enorm));
179:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n", enorm));

181:   /* compare xx and zz. FFTW computes an unnormalized DFT, thus zz = N*x */
182:   s = 1.0 / (PetscReal)N;
183:   PetscCall(VecScale(zz, s));
184:   PetscCall(VecAXPY(xx, -1.0, zz));
185:   PetscCall(VecNorm(xx, NORM_1, &enorm));
186:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n", enorm));

188:   /* compare y and yy: USFFT and FFTW results*/
189:   PetscCall(VecNorm(y, NORM_2, &norm));
190:   PetscCall(VecAXPY(y, -1.0, yy));
191:   PetscCall(VecNorm(y, NORM_1, &enorm));
192:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n", norm));
193:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n", enorm));

195:   /* compare z and zz: USFFT and FFTW results*/
196:   PetscCall(VecNorm(z, NORM_2, &norm));
197:   PetscCall(VecAXPY(z, -1.0, zz));
198:   PetscCall(VecNorm(z, NORM_1, &enorm));
199:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n", norm));
200:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n", enorm));

202:   /* free spaces */
203:   PetscCall(DMRestoreGlobalVector(da, &x));
204:   PetscCall(DMRestoreGlobalVector(da, &xx));
205:   PetscCall(DMRestoreGlobalVector(da, &y));
206:   PetscCall(DMRestoreGlobalVector(da, &yy));
207:   PetscCall(DMRestoreGlobalVector(da, &z));
208:   PetscCall(DMRestoreGlobalVector(da, &zz));

210:   PetscCall(PetscFinalize());
211:   return 0;
212: }