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