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