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 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 28 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex143");CHKERRQ(ierr); 29 ierr = PetscOptionsBool("-vec_view draw", "View the vectors", "ex143", view, &view, NULL);CHKERRQ(ierr); 30 ierr = PetscOptionsBool("-use_FFTW_interface", "Use PETSc-FFTW interface", "ex143",use_interface, &use_interface, NULL);CHKERRQ(ierr); 31 ierr = PetscOptionsEnd();CHKERRQ(ierr); 32 33 ierr = PetscOptionsGetBool(NULL,NULL,"-use_FFTW_interface",&use_interface,NULL);CHKERRQ(ierr); 34 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr); 35 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr); 36 37 ierr = PetscRandomCreate(PETSC_COMM_WORLD, &rdm);CHKERRQ(ierr); 38 ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr); 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) { 49 ierr = PetscPrintf(PETSC_COMM_SELF,"Use FFTW without PETSc-FFTW interface, DIM %D\n",DIM);CHKERRQ(ierr); 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 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,(PetscInt)local_n0*N1,(PetscInt)N,(const PetscScalar*)data_in,&x);CHKERRQ(ierr); 60 ierr = PetscObjectSetName((PetscObject) x, "Real Space vector");CHKERRQ(ierr); 61 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,(PetscInt)local_n0*N1,(PetscInt)N,(const PetscScalar*)data_out,&y);CHKERRQ(ierr); 62 ierr = PetscObjectSetName((PetscObject) y, "Frequency space vector");CHKERRQ(ierr); 63 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,(PetscInt)local_n0*N1,(PetscInt)N,(const PetscScalar*)data_out2,&z);CHKERRQ(ierr); 64 ierr = PetscObjectSetName((PetscObject) z, "Reconstructed vector");CHKERRQ(ierr); 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 ierr = VecSetRandom(x, rdm);CHKERRQ(ierr); 70 if (view) {ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);} 71 72 fftw_execute(fplan); 73 if (view) {ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);} 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 ierr = VecScale(z,a);CHKERRQ(ierr); 80 if (view) {ierr = VecView(z, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);} 81 ierr = VecAXPY(z,-1.0,x);CHKERRQ(ierr); 82 ierr = VecNorm(z,NORM_1,&enorm);CHKERRQ(ierr); 83 if (enorm > 1.e-11 && !rank) { 84 ierr = PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %g\n",(double)enorm);CHKERRQ(ierr); 85 } 86 87 /* Free spaces */ 88 fftw_destroy_plan(fplan); 89 fftw_destroy_plan(bplan); 90 fftw_free(data_in); ierr = VecDestroy(&x);CHKERRQ(ierr); 91 fftw_free(data_out); ierr = VecDestroy(&y);CHKERRQ(ierr); 92 fftw_free(data_out2);ierr = VecDestroy(&z);CHKERRQ(ierr); 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 ierr = PetscMalloc1(i,&dim);CHKERRQ(ierr); 104 for (k=0; k<i; k++) { 105 dim[k]=30; 106 } 107 N *= dim[i-1]; 108 109 110 /* Create FFTW object */ 111 if (!rank) printf("Use PETSc-FFTW interface...%d-DIM: %d\n",(int)DIM,(int)N); 112 113 ierr = MatCreateFFT(PETSC_COMM_WORLD,DIM,dim,MATFFTW,&A);CHKERRQ(ierr); 114 115 /* Create vectors that are compatible with parallel layout of A - must call MatCreateVecs()! */ 116 117 ierr = MatCreateVecsFFTW(A,&x,&y,&z);CHKERRQ(ierr); 118 ierr = PetscObjectSetName((PetscObject) x, "Real space vector");CHKERRQ(ierr); 119 ierr = PetscObjectSetName((PetscObject) y, "Frequency space vector");CHKERRQ(ierr); 120 ierr = PetscObjectSetName((PetscObject) z, "Reconstructed vector");CHKERRQ(ierr); 121 122 /* Set values of space vector x */ 123 ierr = VecSetRandom(x,rdm);CHKERRQ(ierr); 124 125 if (view) {ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);} 126 127 /* Apply FFTW_FORWARD and FFTW_BACKWARD */ 128 ierr = MatMult(A,x,y);CHKERRQ(ierr); 129 if (view) {ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);} 130 131 ierr = MatMultTranspose(A,y,z);CHKERRQ(ierr); 132 133 /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */ 134 a = 1.0/(PetscReal)N; 135 ierr = VecScale(z,a);CHKERRQ(ierr); 136 if (view) {ierr = VecView(z,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);} 137 ierr = VecAXPY(z,-1.0,x);CHKERRQ(ierr); 138 ierr = VecNorm(z,NORM_1,&enorm);CHKERRQ(ierr); 139 if (enorm > 1.e-9 && !rank) { 140 ierr = PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %e\n",enorm);CHKERRQ(ierr); 141 } 142 143 ierr = VecDestroy(&x);CHKERRQ(ierr); 144 ierr = VecDestroy(&y);CHKERRQ(ierr); 145 ierr = VecDestroy(&z);CHKERRQ(ierr); 146 ierr = MatDestroy(&A);CHKERRQ(ierr); 147 148 ierr = PetscFree(dim);CHKERRQ(ierr); 149 } 150 } 151 152 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr); 153 ierr = PetscFinalize(); 154 return ierr; 155 } 156 157 158 /*TEST 159 160 build: 161 requires: fftw complex 162 163 test: 164 output_file: output/ex143.out 165 166 test: 167 suffix: 2 168 nsize: 3 169 output_file: output/ex143.out 170 171 TEST*/ 172