Actual source code: ex146.c
1: /* This program illustrates use of paralllel real FFT*/
2: static char help[]="This program illustrates the use of parallel real 3D fftw (without PETSc interface)";
3: #include <petscmat.h>
4: #include <fftw3.h>
5: #include <fftw3-mpi.h>
7: int main(int argc,char **args)
8: {
9: ptrdiff_t N0=256,N1=256,N2=256,N3=2,dim[4];
10: fftw_plan bplan,fplan;
11: fftw_complex *out;
12: double *in1,*in2;
13: ptrdiff_t alloc_local,local_n0,local_0_start;
14: ptrdiff_t local_n1,local_1_start;
15: PetscInt i,j,indx,n1;
16: PetscInt size,rank,n,N,*in,N_factor,NM;
17: PetscScalar *data_fin,value1,one=1.57,zero=0.0;
18: PetscScalar a,*x_arr,*y_arr,*z_arr,enorm;
19: Vec fin,fout,fout1,ini,final;
20: PetscRandom rnd;
21: VecScatter vecscat,vecscat1;
22: IS indx1,indx2;
23: PetscInt *indx3,k,l,*indx4;
24: PetscInt low,tempindx,tempindx1;
26: PetscInitialize(&argc,&args,(char*)0,help);
27: #if defined(PETSC_USE_COMPLEX)
28: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires real numbers. Your current scalar type is complex");
29: #endif
30: MPI_Comm_size(PETSC_COMM_WORLD, &size);
31: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
33: PetscRandomCreate(PETSC_COMM_WORLD,&rnd);
35: alloc_local = fftw_mpi_local_size_3d_transposed(N0,N1,N2/2+1,PETSC_COMM_WORLD,&local_n0,&local_0_start,&local_n1,&local_1_start);
37: /* printf("The value alloc_local is %ld from process %d\n",alloc_local,rank); */
38: printf("The value local_n0 is %ld from process %d\n",local_n0,rank);
39: /* printf("The value local_0_start is %ld from process %d\n",local_0_start,rank);*/
40: /* printf("The value local_n1 is %ld from process %d\n",local_n1,rank); */
41: /* printf("The value local_1_start is %ld from process %d\n",local_1_start,rank);*/
43: /* Allocate space for input and output arrays */
45: in1=(double*)fftw_malloc(sizeof(double)*alloc_local*2);
46: in2=(double*)fftw_malloc(sizeof(double)*alloc_local*2);
47: out=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);
49: N=2*N0*N1*(N2/2+1);N_factor=N0*N1*N2;
50: n=2*local_n0*N1*(N2/2+1);n1=local_n1*N0*2*N1;
52: /* printf("The value N is %d from process %d\n",N,rank); */
53: /* printf("The value n is %d from process %d\n",n,rank); */
54: /* printf("The value n1 is %d from process %d\n",n1,rank); */
55: /* Creating data vector and accompanying array with VeccreateMPIWithArray */
56: VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)in1,&fin);
57: VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)out,&fout);
58: VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)in2,&fout1);
60: /* VecGetSize(fin,&size); */
61: /* printf("The size is %d\n",size); */
63: VecSet(fin,one);
64: VecSet(fout,zero);
65: VecSet(fout1,zero);
67: VecAssemblyBegin(fin);
68: VecAssemblyEnd(fin);
69: /* VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */
71: VecGetArray(fin,&x_arr);
72: VecGetArray(fout1,&z_arr);
73: VecGetArray(fout,&y_arr);
75: fplan=fftw_mpi_plan_dft_r2c_3d(N0,N1,N2,(double*)x_arr,(fftw_complex*)y_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);
76: bplan=fftw_mpi_plan_dft_c2r_3d(N0,N1,N2,(fftw_complex*)y_arr,(double*)z_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);
78: fftw_execute(fplan);
79: fftw_execute(bplan);
81: VecRestoreArray(fin,&x_arr);
82: VecRestoreArray(fout1,&z_arr);
83: VecRestoreArray(fout,&y_arr);
85: /* a = 1.0/(PetscReal)N_factor; */
86: /* VecScale(fout1,a); */
87: VecCreate(PETSC_COMM_WORLD,&ini);
88: VecCreate(PETSC_COMM_WORLD,&final);
89: VecSetSizes(ini,local_n0*N1*N2,N_factor);
90: VecSetSizes(final,local_n0*N1*N2,N_factor);
91: /* VecSetSizes(ini,PETSC_DECIDE,N_factor); */
92: /* VecSetSizes(final,PETSC_DECIDE,N_factor); */
93: VecSetFromOptions(ini);
94: VecSetFromOptions(final);
96: if (N2%2==0) NM=N2+2;
97: else NM=N2+1;
99: VecGetOwnershipRange(fin,&low,NULL);
100: printf("The local index is %d from %d\n",low,rank);
101: PetscMalloc1(local_n0*N1*N2,&indx3);
102: PetscMalloc1(local_n0*N1*N2,&indx4);
103: for (i=0; i<local_n0; i++) {
104: for (j=0;j<N1;j++) {
105: for (k=0;k<N2;k++) {
106: tempindx = i*N1*N2 + j*N2 + k;
107: tempindx1 = i*N1*NM + j*NM + k;
109: indx3[tempindx]=local_0_start*N1*N2+tempindx;
110: indx4[tempindx]=low+tempindx1;
111: }
112: /* printf("index3 %d from proc %d is \n",indx3[tempindx],rank); */
113: /* printf("index4 %d from proc %d is \n",indx4[tempindx],rank); */
114: }
115: }
116: VecGetValues(fin,local_n0*N1*N2,indx4,x_arr);
117: VecSetValues(ini,local_n0*N1*N2,indx3,x_arr,INSERT_VALUES);
118: VecAssemblyBegin(ini);
119: VecAssemblyEnd(ini);
121: VecGetValues(fout1,local_n0*N1*N2,indx4,y_arr);
122: VecSetValues(final,local_n0*N1*N2,indx3,y_arr,INSERT_VALUES);
123: VecAssemblyBegin(final);
124: VecAssemblyEnd(final);
126: printf("The local index value is %ld from %d",local_n0*N1*N2,rank);
127: /*
128: for (i=0;i<N0;i++) {
129: for (j=0;j<N1;j++) {
130: indx=i*N1*NM+j*NM;
131: ISCreateStride(PETSC_COMM_WORLD,N2,indx,1,&indx1);
132: indx=i*N1*N2+j*N2;
133: ISCreateStride(PETSC_COMM_WORLD,N2,indx,1,&indx2);
134: VecScatterCreate(fin,indx1,ini,indx2,&vecscat);
135: VecScatterBegin(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
136: VecScatterEnd(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
137: VecScatterCreate(fout1,indx1,final,indx2,&vecscat1);
138: VecScatterBegin(vecscat1,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
139: VecScatterEnd(vecscat1,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
140: }
141: }
142: */
143: a = 1.0/(PetscReal)N_factor;
144: VecScale(fout1,a);
145: VecScale(final,a);
147: VecAssemblyBegin(ini);
148: VecAssemblyEnd(ini);
150: VecAssemblyBegin(final);
151: VecAssemblyEnd(final);
153: /* VecView(final,PETSC_VIEWER_STDOUT_WORLD); */
154: VecAXPY(final,-1.0,ini);
155: VecNorm(final,NORM_1,&enorm);
156: PetscPrintf(PETSC_COMM_WORLD," Error norm of |x - z| = %e\n",enorm);
157: fftw_destroy_plan(fplan);
158: fftw_destroy_plan(bplan);
159: fftw_free(in1); VecDestroy(&fin);
160: fftw_free(out); VecDestroy(&fout);
161: fftw_free(in2); VecDestroy(&fout1);
163: PetscFinalize();
164: return 0;
165: }