1 static char help[] = "Illustrate how to use mpi FFTW and PETSc-FFTW interface \n\n"; 2 3 /* 4 Usage: 5 mpiexec -n <np> ./ex158 -use_FFTW_interface NO 6 mpiexec -n <np> ./ex158 -use_FFTW_interface YES 7 */ 8 9 #include <petscmat.h> 10 #include <fftw3-mpi.h> 11 12 int main(int argc,char **args) 13 { 14 PetscMPIInt rank,size; 15 PetscInt N0=50,N1=20,N=N0*N1; 16 PetscRandom rdm; 17 PetscScalar a; 18 PetscReal enorm; 19 Vec x,y,z; 20 PetscBool view=PETSC_FALSE,use_interface=PETSC_TRUE; 21 22 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 23 #if defined(PETSC_USE_COMPLEX) 24 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires real numbers. Your current scalar type is complex"); 25 #endif 26 27 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex158"); 28 PetscCall(PetscOptionsBool("-use_FFTW_interface", "Use PETSc-FFTW interface", "ex158",use_interface, &use_interface, NULL)); 29 PetscOptionsEnd(); 30 31 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 32 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 33 34 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm)); 35 PetscCall(PetscRandomSetFromOptions(rdm)); 36 37 if (!use_interface) { 38 /* Use mpi FFTW without PETSc-FFTW interface, 2D case only */ 39 /*---------------------------------------------------------*/ 40 fftw_plan fplan,bplan; 41 fftw_complex *data_in,*data_out,*data_out2; 42 ptrdiff_t alloc_local,local_n0,local_0_start; 43 44 if (rank == 0) printf("Use FFTW without PETSc-FFTW interface\n"); 45 fftw_mpi_init(); 46 N = N0*N1; 47 alloc_local = fftw_mpi_local_size_2d(N0,N1,PETSC_COMM_WORLD,&local_n0,&local_0_start); 48 49 data_in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 50 data_out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 51 data_out2 = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 52 53 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,(PetscInt)local_n0*N1,(PetscInt)N,(const PetscScalar*)data_in,&x)); 54 PetscCall(PetscObjectSetName((PetscObject) x, "Real Space vector")); 55 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,(PetscInt)local_n0*N1,(PetscInt)N,(const PetscScalar*)data_out,&y)); 56 PetscCall(PetscObjectSetName((PetscObject) y, "Frequency space vector")); 57 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,(PetscInt)local_n0*N1,(PetscInt)N,(const PetscScalar*)data_out2,&z)); 58 PetscCall(PetscObjectSetName((PetscObject) z, "Reconstructed vector")); 59 60 fplan = fftw_mpi_plan_dft_2d(N0,N1,data_in,data_out,PETSC_COMM_WORLD,FFTW_FORWARD,FFTW_ESTIMATE); 61 bplan = fftw_mpi_plan_dft_2d(N0,N1,data_out,data_out2,PETSC_COMM_WORLD,FFTW_BACKWARD,FFTW_ESTIMATE); 62 63 PetscCall(VecSetRandom(x, rdm)); 64 if (view) PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 65 66 fftw_execute(fplan); 67 if (view) PetscCall(VecView(y,PETSC_VIEWER_STDOUT_WORLD)); 68 69 fftw_execute(bplan); 70 71 /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */ 72 a = 1.0/(PetscReal)N; 73 PetscCall(VecScale(z,a)); 74 if (view) PetscCall(VecView(z, PETSC_VIEWER_STDOUT_WORLD)); 75 PetscCall(VecAXPY(z,-1.0,x)); 76 PetscCall(VecNorm(z,NORM_1,&enorm)); 77 if (enorm > 1.e-11) { 78 PetscCall(PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %g\n",(double)enorm)); 79 } 80 81 /* Free spaces */ 82 fftw_destroy_plan(fplan); 83 fftw_destroy_plan(bplan); 84 fftw_free(data_in); PetscCall(VecDestroy(&x)); 85 fftw_free(data_out); PetscCall(VecDestroy(&y)); 86 fftw_free(data_out2);PetscCall(VecDestroy(&z)); 87 88 } else { 89 /* Use PETSc-FFTW interface */ 90 /*-------------------------------------------*/ 91 PetscInt i,*dim,k,DIM; 92 Mat A; 93 Vec input,output; 94 95 N=30; 96 for (i=2; i<3; i++) { /* (i=3,4: -- error in VecScatterPetscToFFTW(A,input,x); */ 97 DIM = i; 98 PetscCall(PetscMalloc1(i,&dim)); 99 for (k=0; k<i; k++) { 100 dim[k]=30; 101 } 102 N *= dim[i-1]; 103 104 /* Create FFTW object */ 105 if (rank == 0) { 106 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Use PETSc-FFTW interface...%d-DIM:%d \n",DIM,N)); 107 } 108 PetscCall(MatCreateFFT(PETSC_COMM_WORLD,DIM,dim,MATFFTW,&A)); 109 110 /* Create FFTW vectors that are compatible with parallel layout of A */ 111 PetscCall(MatCreateVecsFFTW(A,&x,&y,&z)); 112 PetscCall(PetscObjectSetName((PetscObject) x, "Real space vector")); 113 PetscCall(PetscObjectSetName((PetscObject) y, "Frequency space vector")); 114 PetscCall(PetscObjectSetName((PetscObject) z, "Reconstructed vector")); 115 116 /* Create and set PETSc vector */ 117 PetscCall(VecCreate(PETSC_COMM_WORLD,&input)); 118 PetscCall(VecSetSizes(input,PETSC_DECIDE,N)); 119 PetscCall(VecSetFromOptions(input)); 120 PetscCall(VecSetRandom(input,rdm)); 121 PetscCall(VecDuplicate(input,&output)); 122 if (view) PetscCall(VecView(input,PETSC_VIEWER_STDOUT_WORLD)); 123 124 /* Vector input is copied to another vector x using VecScatterPetscToFFTW. This is because the user data 125 can have any parallel layout. But FFTW requires special parallel layout of the data. Hence the original 126 data which is in the vector "input" here, needs to be copied to a vector x, which has the correct parallel 127 layout for FFTW. Also, during parallel real transform, this pads extra zeros automatically 128 at the end of last dimension. This padding is required by FFTW to perform parallel real D.F.T. */ 129 PetscCall(VecScatterPetscToFFTW(A,input,x));/* buggy for dim = 3, 4... */ 130 131 /* Apply FFTW_FORWARD and FFTW_BACKWARD */ 132 PetscCall(MatMult(A,x,y)); 133 if (view) PetscCall(VecView(y,PETSC_VIEWER_STDOUT_WORLD)); 134 PetscCall(MatMultTranspose(A,y,z)); 135 136 /* Output from Backward DFT needs to be modified to obtain user readable data the routine VecScatterFFTWToPetsc 137 performs the job. In some sense this is the reverse operation of VecScatterPetscToFFTW. This routine gets rid of 138 the extra spaces that were artificially padded to perform real parallel transform. */ 139 PetscCall(VecScatterFFTWToPetsc(A,z,output)); 140 141 /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */ 142 a = 1.0/(PetscReal)N; 143 PetscCall(VecScale(output,a)); 144 if (view) PetscCall(VecView(output,PETSC_VIEWER_STDOUT_WORLD)); 145 PetscCall(VecAXPY(output,-1.0,input)); 146 PetscCall(VecNorm(output,NORM_1,&enorm)); 147 if (enorm > 1.e-09 && rank == 0) { 148 PetscCall(PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %e\n",enorm)); 149 } 150 151 /* Free spaces */ 152 PetscCall(PetscFree(dim)); 153 PetscCall(VecDestroy(&input)); 154 PetscCall(VecDestroy(&output)); 155 PetscCall(VecDestroy(&x)); 156 PetscCall(VecDestroy(&y)); 157 PetscCall(VecDestroy(&z)); 158 PetscCall(MatDestroy(&A)); 159 } 160 } 161 PetscCall(PetscRandomDestroy(&rdm)); 162 PetscCall(PetscFinalize()); 163 return 0; 164 } 165 166 /*TEST 167 168 build: 169 requires: !mpiuni fftw !complex 170 171 test: 172 output_file: output/ex158.out 173 174 test: 175 suffix: 2 176 nsize: 3 177 178 TEST*/ 179