Actual source code: ex27.c
petsc-3.9.4 2018-09-11
1: static char help[] = "Test sequential USFFT interface on a uniform DMDA and compares the result to FFTW\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: #include <petscmat.h>
11: #include <petscdm.h>
12: #include <petscdmda.h>
13: PetscInt main(PetscInt 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;
32: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
33: MPI_Comm_size(PETSC_COMM_WORLD, &size);
34: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This is a uniprocessor example only!");
35: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "USFFT Options", "ex27");
36: PetscOptionsEList("-function", "Function type", "ex27", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
37: function = (FuncType) func;
38: PetscOptionsEnd();
39: PetscOptionsGetBool(NULL,NULL,"-view_x",&view_x,NULL);
40: PetscOptionsGetBool(NULL,NULL,"-view_y",&view_y,NULL);
41: PetscOptionsGetBool(NULL,NULL,"-view_z",&view_z,NULL);
42: PetscOptionsGetIntArray(NULL,NULL,"-dim",dim,&ndim,NULL);
44: 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: DMSetFromOptions(da);
47: DMSetUp(da);
49: /* Coordinates */
50: DMGetCoordinateDM(da, &coordsda);
51: DMGetGlobalVector(coordsda, &coords);
52: 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: 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: VecRestoreArray(coords, &a);
68: }
69: DMSetCoordinates(da, coords);
71: /* Work vectors */
72: DMGetGlobalVector(da, &x);
73: PetscObjectSetName((PetscObject) x, "Real space vector");
74: DMGetGlobalVector(da, &xx);
75: PetscObjectSetName((PetscObject) xx, "Real space vector");
76: DMGetGlobalVector(da, &y);
77: PetscObjectSetName((PetscObject) y, "USFFT frequency space vector");
78: DMGetGlobalVector(da, &yy);
79: PetscObjectSetName((PetscObject) yy, "FFTW frequency space vector");
80: DMGetGlobalVector(da, &z);
81: PetscObjectSetName((PetscObject) z, "USFFT reconstructed vector");
82: DMGetGlobalVector(da, &zz);
83: PetscObjectSetName((PetscObject) zz, "FFTW reconstructed vector");
85: PetscPrintf(PETSC_COMM_SELF, "%3-D: USFFT on vector of ");
86: for (i = 0, N = 1; i < 3; i++) {
87: PetscPrintf(PETSC_COMM_SELF, "dim[%d] = %d ",i,dim[i]);
88: N *= dim[i];
89: }
90: PetscPrintf(PETSC_COMM_SELF, "; total size %d \n",N);
93: if (function == RANDOM) {
94: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
95: PetscRandomSetFromOptions(rdm);
96: VecSetRandom(x, rdm);
97: PetscRandomDestroy(&rdm);
98: } else if (function == CONSTANT) {
99: VecSet(x, 1.0);
100: } else if (function == TANH) {
101: PetscScalar *a;
102: VecGetArray(x, &a);
103: PetscInt j,k = 0;
104: for (i = 0; i < 3; ++i) {
105: for (j = 0; j < dim[i]; ++j) {
106: a[k] = tanh((j - dim[i]/2.0)*(10.0/dim[i]));
107: ++k;
108: }
109: }
110: VecRestoreArray(x, &a);
111: }
112: if (view_x) {
113: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
114: }
115: VecCopy(x,xx);
117: VecNorm(x,NORM_2,&norm);
118: PetscPrintf(PETSC_COMM_SELF, "|x|_2 = %g\n",norm);
120: /* create USFFT object */
121: MatCreateSeqUSFFT(coords,da,&A);
122: /* create FFTW object */
123: MatCreateSeqFFTW(PETSC_COMM_SELF,3,dim,&AA);
125: /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */
126: MatMult(A,x,z);
127: MatMult(AA,xx,zz);
128: /* Now apply USFFT and FFTW forward several (3) times */
129: for (i=0; i<3; ++i) {
130: MatMult(A,x,y);
131: MatMult(AA,xx,yy);
132: MatMultTranspose(A,y,z);
133: MatMultTranspose(AA,yy,zz);
134: }
136: if (view_y) {
137: PetscPrintf(PETSC_COMM_WORLD, "y = \n");
138: VecView(y, PETSC_VIEWER_STDOUT_WORLD);
139: PetscPrintf(PETSC_COMM_WORLD, "yy = \n");
140: VecView(yy, PETSC_VIEWER_STDOUT_WORLD);
141: }
143: if (view_z) {
144: PetscPrintf(PETSC_COMM_WORLD, "z = \n");
145: VecView(z, PETSC_VIEWER_STDOUT_WORLD);
146: PetscPrintf(PETSC_COMM_WORLD, "zz = \n");
147: VecView(zz, PETSC_VIEWER_STDOUT_WORLD);
148: }
150: /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */
151: s = 1.0/(PetscReal)N;
152: VecScale(z,s);
153: VecAXPY(x,-1.0,z);
154: VecNorm(x,NORM_1,&enorm);
155: PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n",enorm);
157: /* compare xx and zz. FFTW computes an unnormalized DFT, thus zz = N*x */
158: s = 1.0/(PetscReal)N;
159: VecScale(zz,s);
160: VecAXPY(xx,-1.0,zz);
161: VecNorm(xx,NORM_1,&enorm);
162: PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n",enorm);
164: /* compare y and yy: USFFT and FFTW results*/
165: VecNorm(y,NORM_2,&norm);
166: VecAXPY(y,-1.0,yy);
167: VecNorm(y,NORM_1,&enorm);
168: PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n",norm);
169: PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n",enorm);
171: /* compare z and zz: USFFT and FFTW results*/
172: VecNorm(z,NORM_2,&norm);
173: VecAXPY(z,-1.0,zz);
174: VecNorm(z,NORM_1,&enorm);
175: PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n",norm);
176: PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n",enorm);
179: /* free spaces */
180: DMRestoreGlobalVector(da,&x);
181: DMRestoreGlobalVector(da,&xx);
182: DMRestoreGlobalVector(da,&y);
183: DMRestoreGlobalVector(da,&yy);
184: DMRestoreGlobalVector(da,&z);
185: DMRestoreGlobalVector(da,&zz);
186: VecDestroy(&coords);
187: DMDestroy(&da);
188: PetscFinalize();
189: return ierr;
190: }