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>
main(PetscInt argc,char ** args)15 PetscInt main(PetscInt 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);CHKERRQ(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);
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 
105   ierr = PetscPrintf(PETSC_COMM_SELF, "%3-D: USFFT on vector of ");CHKERRQ(ierr);
106   for (i = 0, N = 1; i < 3; i++) {
107     ierr = PetscPrintf(PETSC_COMM_SELF, "dim[%d] = %d ",i,dim[i]);CHKERRQ(ierr);
108     N   *= dim[i];
109   }
110   ierr = PetscPrintf(PETSC_COMM_SELF, "; total size %d \n",N);CHKERRQ(ierr);
111 
112 
113   if (function == RANDOM) {
114     ierr = PetscRandomCreate(PETSC_COMM_SELF, &rdm);CHKERRQ(ierr);
115     ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr);
116     ierr = VecSetRandom(x, rdm);CHKERRQ(ierr);
117     ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr);
118   } else if (function == CONSTANT) {
119     ierr = VecSet(x, 1.0);CHKERRQ(ierr);
120   } else if (function == TANH) {
121     PetscScalar *a;
122     ierr = VecGetArray(x, &a);CHKERRQ(ierr);
123     PetscInt j,k = 0;
124     for (i = 0; i < 3; ++i) {
125       for (j = 0; j < dim[i]; ++j) {
126         a[k] = tanh((j - dim[i]/2.0)*(10.0/dim[i]));
127         ++k;
128       }
129     }
130     ierr = VecRestoreArray(x, &a);CHKERRQ(ierr);
131   }
132   if (view_x) {
133     ierr = VecView(x, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
134   }
135   ierr = VecCopy(x,xx);CHKERRQ(ierr);
136   /* Split xx */
137   ierr = VecStrideGatherAll(xx,xxsplit, INSERT_VALUES);CHKERRQ(ierr); /*YES! 'Gather' means 'split' (or maybe 'scatter'?)! */
138 
139   ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
140   ierr = PetscPrintf(PETSC_COMM_SELF, "|x|_2 = %g\n",norm);CHKERRQ(ierr);
141 
142   /* create USFFT object */
143   ierr = MatCreateSeqUSFFT(da,da,&A);CHKERRQ(ierr);
144   /* create FFTW object */
145   ierr = MatCreateSeqFFTW(PETSC_COMM_SELF,3,dim,&AA);CHKERRQ(ierr);
146 
147   /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */
148   ierr = MatMult(A,x,z);CHKERRQ(ierr);
149   for (ii = 0; ii < 3; ++ii) {
150     ierr = MatMult(AA,xxsplit[ii],zzsplit[ii]);CHKERRQ(ierr);
151   }
152   /* Now apply USFFT and FFTW forward several (3) times */
153   for (i=0; i<3; ++i) {
154     ierr = MatMult(A,x,y);CHKERRQ(ierr);
155     for (ii = 0; ii < 3; ++ii) {
156       ierr = MatMult(AA,xxsplit[ii],yysplit[ii]);CHKERRQ(ierr);
157     }
158     ierr = MatMultTranspose(A,y,z);CHKERRQ(ierr);
159     for (ii = 0; ii < 3; ++ii) {
160       ierr = MatMult(AA,yysplit[ii],zzsplit[ii]);CHKERRQ(ierr);
161     }
162   }
163   /* Unsplit yy */
164   ierr = VecStrideScatterAll(yysplit, yy, INSERT_VALUES);CHKERRQ(ierr); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */
165   /* Unsplit zz */
166   ierr = VecStrideScatterAll(zzsplit, zz, INSERT_VALUES);CHKERRQ(ierr); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */
167 
168   if (view_y) {
169     ierr = PetscPrintf(PETSC_COMM_WORLD, "y = \n");CHKERRQ(ierr);
170     ierr = VecView(y, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
171     ierr = PetscPrintf(PETSC_COMM_WORLD, "yy = \n");CHKERRQ(ierr);
172     ierr = VecView(yy, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
173   }
174 
175   if (view_z) {
176     ierr = PetscPrintf(PETSC_COMM_WORLD, "z = \n");CHKERRQ(ierr);
177     ierr = VecView(z, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
178     ierr = PetscPrintf(PETSC_COMM_WORLD, "zz = \n");CHKERRQ(ierr);
179     ierr = VecView(zz, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
180   }
181 
182   /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */
183   s    = 1.0/(PetscReal)N;
184   ierr = VecScale(z,s);CHKERRQ(ierr);
185   ierr = VecAXPY(x,-1.0,z);CHKERRQ(ierr);
186   ierr = VecNorm(x,NORM_1,&enorm);CHKERRQ(ierr);
187   ierr = PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n",enorm);CHKERRQ(ierr);
188 
189   /* compare xx and zz. FFTW computes an unnormalized DFT, thus zz = N*x */
190   s    = 1.0/(PetscReal)N;
191   ierr = VecScale(zz,s);CHKERRQ(ierr);
192   ierr = VecAXPY(xx,-1.0,zz);CHKERRQ(ierr);
193   ierr = VecNorm(xx,NORM_1,&enorm);CHKERRQ(ierr);
194   ierr = PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n",enorm);CHKERRQ(ierr);
195 
196   /* compare y and yy: USFFT and FFTW results*/
197   ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr);
198   ierr = VecAXPY(y,-1.0,yy);CHKERRQ(ierr);
199   ierr = VecNorm(y,NORM_1,&enorm);CHKERRQ(ierr);
200   ierr = PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n",norm);CHKERRQ(ierr);
201   ierr = PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n",enorm);CHKERRQ(ierr);
202 
203   /* compare z and zz: USFFT and FFTW results*/
204   ierr = VecNorm(z,NORM_2,&norm);CHKERRQ(ierr);
205   ierr = VecAXPY(z,-1.0,zz);CHKERRQ(ierr);
206   ierr = VecNorm(z,NORM_1,&enorm);CHKERRQ(ierr);
207   ierr = PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n",norm);CHKERRQ(ierr);
208   ierr = PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n",enorm);CHKERRQ(ierr);
209 
210 
211   /* free spaces */
212   ierr = DMRestoreGlobalVector(da,&x);CHKERRQ(ierr);
213   ierr = DMRestoreGlobalVector(da,&xx);CHKERRQ(ierr);
214   ierr = DMRestoreGlobalVector(da,&y);CHKERRQ(ierr);
215   ierr = DMRestoreGlobalVector(da,&yy);CHKERRQ(ierr);
216   ierr = DMRestoreGlobalVector(da,&z);CHKERRQ(ierr);
217   ierr = DMRestoreGlobalVector(da,&zz);CHKERRQ(ierr);
218 
219   ierr = PetscFinalize();
220   return ierr;
221 }
222