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