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