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