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