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"; 2 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. 7 8 */ 9 10 #define DOF 3 11 12 #include <petscmat.h> 13 #include <petscdm.h> 14 #include <petscdmda.h> 15 int main(int argc, char **args) { 16 typedef enum { 17 RANDOM, 18 CONSTANT, 19 TANH, 20 NUM_FUNCS 21 } FuncType; 22 const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"}; 23 Mat A, AA; 24 PetscMPIInt size; 25 PetscInt N, i, stencil = 1, dof = 3; 26 PetscInt dim[3] = {10, 10, 10}, ndim = 3; 27 Vec coords, x, y, z, xx, yy, zz; 28 Vec xxsplit[DOF], yysplit[DOF], zzsplit[DOF]; 29 PetscReal h[3]; 30 PetscScalar s; 31 PetscRandom rdm; 32 PetscReal norm, enorm; 33 PetscInt func, ii; 34 FuncType function = TANH; 35 DM da, da1, coordsda; 36 PetscBool view_x = PETSC_FALSE, view_y = PETSC_FALSE, view_z = PETSC_FALSE; 37 38 PetscFunctionBeginUser; 39 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 40 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 41 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only!"); 42 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "USFFT Options", "ex27"); 43 PetscCall(PetscOptionsEList("-function", "Function type", "ex27", funcNames, NUM_FUNCS, funcNames[function], &func, NULL)); 44 function = (FuncType)func; 45 PetscOptionsEnd(); 46 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_x", &view_x, NULL)); 47 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_y", &view_y, NULL)); 48 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_z", &view_z, NULL)); 49 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dim", dim, &ndim, NULL)); 50 51 /* DMDA with the correct fiber dimension */ 52 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)); 53 PetscCall(DMSetFromOptions(da)); 54 PetscCall(DMSetUp(da)); 55 /* DMDA with fiber dimension 1 for split fields */ 56 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)); 57 PetscCall(DMSetFromOptions(da1)); 58 PetscCall(DMSetUp(da1)); 59 60 /* Coordinates */ 61 PetscCall(DMGetCoordinateDM(da, &coordsda)); 62 PetscCall(DMGetGlobalVector(coordsda, &coords)); 63 PetscCall(PetscObjectSetName((PetscObject)coords, "Grid coordinates")); 64 for (i = 0, N = 1; i < 3; i++) { 65 h[i] = 1.0 / dim[i]; 66 PetscScalar *a; 67 PetscCall(VecGetArray(coords, &a)); 68 PetscInt j, k, n = 0; 69 for (i = 0; i < 3; ++i) { 70 for (j = 0; j < dim[i]; ++j) { 71 for (k = 0; k < 3; ++k) { 72 a[n] = j * h[i]; /* coordinate along the j-th point in the i-th dimension */ 73 ++n; 74 } 75 } 76 } 77 PetscCall(VecRestoreArray(coords, &a)); 78 } 79 PetscCall(DMSetCoordinates(da, coords)); 80 PetscCall(VecDestroy(&coords)); 81 82 /* Work vectors */ 83 PetscCall(DMGetGlobalVector(da, &x)); 84 PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector")); 85 PetscCall(DMGetGlobalVector(da, &xx)); 86 PetscCall(PetscObjectSetName((PetscObject)xx, "Real space vector")); 87 PetscCall(DMGetGlobalVector(da, &y)); 88 PetscCall(PetscObjectSetName((PetscObject)y, "USFFT frequency space vector")); 89 PetscCall(DMGetGlobalVector(da, &yy)); 90 PetscCall(PetscObjectSetName((PetscObject)yy, "FFTW frequency space vector")); 91 PetscCall(DMGetGlobalVector(da, &z)); 92 PetscCall(PetscObjectSetName((PetscObject)z, "USFFT reconstructed vector")); 93 PetscCall(DMGetGlobalVector(da, &zz)); 94 PetscCall(PetscObjectSetName((PetscObject)zz, "FFTW reconstructed vector")); 95 /* Split vectors for FFTW */ 96 for (ii = 0; ii < 3; ++ii) { 97 PetscCall(DMGetGlobalVector(da1, &xxsplit[ii])); 98 PetscCall(PetscObjectSetName((PetscObject)xxsplit[ii], "Real space split vector")); 99 PetscCall(DMGetGlobalVector(da1, &yysplit[ii])); 100 PetscCall(PetscObjectSetName((PetscObject)yysplit[ii], "FFTW frequency space split vector")); 101 PetscCall(DMGetGlobalVector(da1, &zzsplit[ii])); 102 PetscCall(PetscObjectSetName((PetscObject)zzsplit[ii], "FFTW reconstructed split vector")); 103 } 104 105 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%3-" PetscInt_FMT ": USFFT on vector of ")); 106 for (i = 0, N = 1; i < 3; i++) { 107 PetscCall(PetscPrintf(PETSC_COMM_SELF, "dim[%d] = %d ", i, dim[i])); 108 N *= dim[i]; 109 } 110 PetscCall(PetscPrintf(PETSC_COMM_SELF, "; total size %d \n", N)); 111 112 if (function == RANDOM) { 113 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm)); 114 PetscCall(PetscRandomSetFromOptions(rdm)); 115 PetscCall(VecSetRandom(x, rdm)); 116 PetscCall(PetscRandomDestroy(&rdm)); 117 } else if (function == CONSTANT) { 118 PetscCall(VecSet(x, 1.0)); 119 } else if (function == TANH) { 120 PetscScalar *a; 121 PetscCall(VecGetArray(x, &a)); 122 PetscInt j, k = 0; 123 for (i = 0; i < 3; ++i) { 124 for (j = 0; j < dim[i]; ++j) { 125 a[k] = tanh((j - dim[i] / 2.0) * (10.0 / dim[i])); 126 ++k; 127 } 128 } 129 PetscCall(VecRestoreArray(x, &a)); 130 } 131 if (view_x) PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 132 PetscCall(VecCopy(x, xx)); 133 /* Split xx */ 134 PetscCall(VecStrideGatherAll(xx, xxsplit, INSERT_VALUES)); /*YES! 'Gather' means 'split' (or maybe 'scatter'?)! */ 135 136 PetscCall(VecNorm(x, NORM_2, &norm)); 137 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|x|_2 = %g\n", norm)); 138 139 /* create USFFT object */ 140 PetscCall(MatCreateSeqUSFFT(da, da, &A)); 141 /* create FFTW object */ 142 PetscCall(MatCreateSeqFFTW(PETSC_COMM_SELF, 3, dim, &AA)); 143 144 /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */ 145 PetscCall(MatMult(A, x, z)); 146 for (ii = 0; ii < 3; ++ii) { PetscCall(MatMult(AA, xxsplit[ii], zzsplit[ii])); } 147 /* Now apply USFFT and FFTW forward several (3) times */ 148 for (i = 0; i < 3; ++i) { 149 PetscCall(MatMult(A, x, y)); 150 for (ii = 0; ii < 3; ++ii) { PetscCall(MatMult(AA, xxsplit[ii], yysplit[ii])); } 151 PetscCall(MatMultTranspose(A, y, z)); 152 for (ii = 0; ii < 3; ++ii) { PetscCall(MatMult(AA, yysplit[ii], zzsplit[ii])); } 153 } 154 /* Unsplit yy */ 155 PetscCall(VecStrideScatterAll(yysplit, yy, INSERT_VALUES)); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */ 156 /* Unsplit zz */ 157 PetscCall(VecStrideScatterAll(zzsplit, zz, INSERT_VALUES)); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */ 158 159 if (view_y) { 160 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "y = \n")); 161 PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD)); 162 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "yy = \n")); 163 PetscCall(VecView(yy, PETSC_VIEWER_STDOUT_WORLD)); 164 } 165 166 if (view_z) { 167 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "z = \n")); 168 PetscCall(VecView(z, PETSC_VIEWER_STDOUT_WORLD)); 169 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "zz = \n")); 170 PetscCall(VecView(zz, PETSC_VIEWER_STDOUT_WORLD)); 171 } 172 173 /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */ 174 s = 1.0 / (PetscReal)N; 175 PetscCall(VecScale(z, s)); 176 PetscCall(VecAXPY(x, -1.0, z)); 177 PetscCall(VecNorm(x, NORM_1, &enorm)); 178 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n", enorm)); 179 180 /* compare xx and zz. FFTW computes an unnormalized DFT, thus zz = N*x */ 181 s = 1.0 / (PetscReal)N; 182 PetscCall(VecScale(zz, s)); 183 PetscCall(VecAXPY(xx, -1.0, zz)); 184 PetscCall(VecNorm(xx, NORM_1, &enorm)); 185 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n", enorm)); 186 187 /* compare y and yy: USFFT and FFTW results*/ 188 PetscCall(VecNorm(y, NORM_2, &norm)); 189 PetscCall(VecAXPY(y, -1.0, yy)); 190 PetscCall(VecNorm(y, NORM_1, &enorm)); 191 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n", norm)); 192 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n", enorm)); 193 194 /* compare z and zz: USFFT and FFTW results*/ 195 PetscCall(VecNorm(z, NORM_2, &norm)); 196 PetscCall(VecAXPY(z, -1.0, zz)); 197 PetscCall(VecNorm(z, NORM_1, &enorm)); 198 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n", norm)); 199 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n", enorm)); 200 201 /* free spaces */ 202 PetscCall(DMRestoreGlobalVector(da, &x)); 203 PetscCall(DMRestoreGlobalVector(da, &xx)); 204 PetscCall(DMRestoreGlobalVector(da, &y)); 205 PetscCall(DMRestoreGlobalVector(da, &yy)); 206 PetscCall(DMRestoreGlobalVector(da, &z)); 207 PetscCall(DMRestoreGlobalVector(da, &zz)); 208 209 PetscCall(PetscFinalize()); 210 return 0; 211 } 212