1 static char help[] = "Illustrate how to use mpi FFTW and PETSc-FFTW interface \n\n"; 2 3 /* 4 Compiling the code: 5 This code uses the complex numbers version of PETSc, so configure 6 must be run to enable this 7 8 Usage: 9 mpiexec -n <np> ./ex143 -use_FFTW_interface NO 10 mpiexec -n <np> ./ex143 -use_FFTW_interface YES 11 */ 12 13 #include <petscmat.h> 14 #include <fftw3-mpi.h> 15 16 int main(int argc,char **args) 17 { 18 PetscMPIInt rank,size; 19 PetscInt N0=50,N1=20,N=N0*N1,DIM; 20 PetscRandom rdm; 21 PetscScalar a; 22 PetscReal enorm; 23 Vec x,y,z; 24 PetscBool view=PETSC_FALSE,use_interface=PETSC_TRUE; 25 26 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 27 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex143"); 28 PetscCall(PetscOptionsBool("-vec_view draw", "View the vectors", "ex143", view, &view, NULL)); 29 PetscCall(PetscOptionsBool("-use_FFTW_interface", "Use PETSc-FFTW interface", "ex143",use_interface, &use_interface, NULL)); 30 PetscOptionsEnd(); 31 32 PetscCall(PetscOptionsGetBool(NULL,NULL,"-use_FFTW_interface",&use_interface,NULL)); 33 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 34 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 35 36 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm)); 37 PetscCall(PetscRandomSetFromOptions(rdm)); 38 39 if (!use_interface) { 40 /* Use mpi FFTW without PETSc-FFTW interface, 2D case only */ 41 /*---------------------------------------------------------*/ 42 fftw_plan fplan,bplan; 43 fftw_complex *data_in,*data_out,*data_out2; 44 ptrdiff_t alloc_local,local_n0,local_0_start; 45 46 DIM = 2; 47 if (rank == 0) { 48 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Use FFTW without PETSc-FFTW interface, DIM %" PetscInt_FMT "\n",DIM)); 49 } 50 fftw_mpi_init(); 51 N = N0*N1; 52 alloc_local = fftw_mpi_local_size_2d(N0,N1,PETSC_COMM_WORLD,&local_n0,&local_0_start); 53 54 data_in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 55 data_out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 56 data_out2 = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 57 58 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,(PetscInt)local_n0*N1,(PetscInt)N,(const PetscScalar*)data_in,&x)); 59 PetscCall(PetscObjectSetName((PetscObject) x, "Real Space vector")); 60 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,(PetscInt)local_n0*N1,(PetscInt)N,(const PetscScalar*)data_out,&y)); 61 PetscCall(PetscObjectSetName((PetscObject) y, "Frequency space vector")); 62 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,(PetscInt)local_n0*N1,(PetscInt)N,(const PetscScalar*)data_out2,&z)); 63 PetscCall(PetscObjectSetName((PetscObject) z, "Reconstructed vector")); 64 65 fplan = fftw_mpi_plan_dft_2d(N0,N1,data_in,data_out,PETSC_COMM_WORLD,FFTW_FORWARD,FFTW_ESTIMATE); 66 bplan = fftw_mpi_plan_dft_2d(N0,N1,data_out,data_out2,PETSC_COMM_WORLD,FFTW_BACKWARD,FFTW_ESTIMATE); 67 68 PetscCall(VecSetRandom(x, rdm)); 69 if (view) PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 70 71 fftw_execute(fplan); 72 if (view) PetscCall(VecView(y,PETSC_VIEWER_STDOUT_WORLD)); 73 74 fftw_execute(bplan); 75 76 /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */ 77 a = 1.0/(PetscReal)N; 78 PetscCall(VecScale(z,a)); 79 if (view) PetscCall(VecView(z, PETSC_VIEWER_STDOUT_WORLD)); 80 PetscCall(VecAXPY(z,-1.0,x)); 81 PetscCall(VecNorm(z,NORM_1,&enorm)); 82 if (enorm > 1.e-11 && rank == 0) { 83 PetscCall(PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %g\n",(double)enorm)); 84 } 85 86 /* Free spaces */ 87 fftw_destroy_plan(fplan); 88 fftw_destroy_plan(bplan); 89 fftw_free(data_in); PetscCall(VecDestroy(&x)); 90 fftw_free(data_out); PetscCall(VecDestroy(&y)); 91 fftw_free(data_out2);PetscCall(VecDestroy(&z)); 92 93 } else { 94 /* Use PETSc-FFTW interface */ 95 /*-------------------------------------------*/ 96 PetscInt i,*dim,k; 97 Mat A; 98 99 N=1; 100 for (i=1; i<5; i++) { 101 DIM = i; 102 PetscCall(PetscMalloc1(i,&dim)); 103 for (k=0; k<i; k++) { 104 dim[k]=30; 105 } 106 N *= dim[i-1]; 107 108 /* Create FFTW object */ 109 if (rank == 0) printf("Use PETSc-FFTW interface...%d-DIM: %d\n",(int)DIM,(int)N); 110 111 PetscCall(MatCreateFFT(PETSC_COMM_WORLD,DIM,dim,MATFFTW,&A)); 112 113 /* Create vectors that are compatible with parallel layout of A - must call MatCreateVecs()! */ 114 115 PetscCall(MatCreateVecsFFTW(A,&x,&y,&z)); 116 PetscCall(PetscObjectSetName((PetscObject) x, "Real space vector")); 117 PetscCall(PetscObjectSetName((PetscObject) y, "Frequency space vector")); 118 PetscCall(PetscObjectSetName((PetscObject) z, "Reconstructed vector")); 119 120 /* Set values of space vector x */ 121 PetscCall(VecSetRandom(x,rdm)); 122 123 if (view) PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 124 125 /* Apply FFTW_FORWARD and FFTW_BACKWARD */ 126 PetscCall(MatMult(A,x,y)); 127 if (view) PetscCall(VecView(y,PETSC_VIEWER_STDOUT_WORLD)); 128 129 PetscCall(MatMultTranspose(A,y,z)); 130 131 /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */ 132 a = 1.0/(PetscReal)N; 133 PetscCall(VecScale(z,a)); 134 if (view) PetscCall(VecView(z,PETSC_VIEWER_STDOUT_WORLD)); 135 PetscCall(VecAXPY(z,-1.0,x)); 136 PetscCall(VecNorm(z,NORM_1,&enorm)); 137 if (enorm > 1.e-9 && rank == 0) { 138 PetscCall(PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %e\n",enorm)); 139 } 140 141 PetscCall(VecDestroy(&x)); 142 PetscCall(VecDestroy(&y)); 143 PetscCall(VecDestroy(&z)); 144 PetscCall(MatDestroy(&A)); 145 146 PetscCall(PetscFree(dim)); 147 } 148 } 149 150 PetscCall(PetscRandomDestroy(&rdm)); 151 PetscCall(PetscFinalize()); 152 return 0; 153 } 154 155 /*TEST 156 157 build: 158 requires: !mpiuni fftw complex 159 160 test: 161 output_file: output/ex143.out 162 163 test: 164 suffix: 2 165 nsize: 3 166 output_file: output/ex143.out 167 168 TEST*/ 169