1 static char help[] = "Test sequential USFFT interface on a uniform DMDA and compares the result to FFTW\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 #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; 35 36 PetscFunctionBeginUser; 37 PetscCall(PetscInitialize(&argc, &args, NULL, help)); 38 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 39 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "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)); 48 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)); 52 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)); 73 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")); 87 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)); 94 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)); 116 117 PetscCall(VecNorm(x, NORM_2, &norm)); 118 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|x|_2 = %g\n", norm)); 119 120 /* create USFFT object */ 121 PetscCall(MatCreateSeqUSFFT(coords, da, &A)); 122 /* create FFTW object */ 123 PetscCall(MatCreateSeqFFTW(PETSC_COMM_SELF, 3, dim, &AA)); 124 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 } 135 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 } 142 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 } 149 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)); 156 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)); 163 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)); 170 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)); 177 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 } 190