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 PetscFunctionBeginUser; 27 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 28 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex143"); 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 PetscOptionsEnd(); 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) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Use FFTW without PETSc-FFTW interface, DIM %" PetscInt_FMT "\n", DIM)); 49 fftw_mpi_init(); 50 N = N0 * N1; 51 alloc_local = fftw_mpi_local_size_2d(N0, N1, PETSC_COMM_WORLD, &local_n0, &local_0_start); 52 53 data_in = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * alloc_local); 54 data_out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * alloc_local); 55 data_out2 = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * alloc_local); 56 57 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, (PetscInt)local_n0 * N1, (PetscInt)N, (const PetscScalar *)data_in, &x)); 58 PetscCall(PetscObjectSetName((PetscObject)x, "Real Space vector")); 59 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, (PetscInt)local_n0 * N1, (PetscInt)N, (const PetscScalar *)data_out, &y)); 60 PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector")); 61 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, (PetscInt)local_n0 * N1, (PetscInt)N, (const PetscScalar *)data_out2, &z)); 62 PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector")); 63 64 fplan = fftw_mpi_plan_dft_2d(N0, N1, data_in, data_out, PETSC_COMM_WORLD, FFTW_FORWARD, FFTW_ESTIMATE); 65 bplan = fftw_mpi_plan_dft_2d(N0, N1, data_out, data_out2, PETSC_COMM_WORLD, FFTW_BACKWARD, FFTW_ESTIMATE); 66 67 PetscCall(VecSetRandom(x, rdm)); 68 if (view) PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 69 70 fftw_execute(fplan); 71 if (view) PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD)); 72 73 fftw_execute(bplan); 74 75 /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */ 76 a = 1.0 / (PetscReal)N; 77 PetscCall(VecScale(z, a)); 78 if (view) PetscCall(VecView(z, PETSC_VIEWER_STDOUT_WORLD)); 79 PetscCall(VecAXPY(z, -1.0, x)); 80 PetscCall(VecNorm(z, NORM_1, &enorm)); 81 if (enorm > 1.e-11 && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %g\n", (double)enorm)); 82 83 /* Free spaces */ 84 fftw_destroy_plan(fplan); 85 fftw_destroy_plan(bplan); 86 fftw_free(data_in); 87 PetscCall(VecDestroy(&x)); 88 fftw_free(data_out); 89 PetscCall(VecDestroy(&y)); 90 fftw_free(data_out2); 91 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++) dim[k] = 30; 104 N *= dim[i - 1]; 105 106 /* Create FFTW object */ 107 if (rank == 0) printf("Use PETSc-FFTW interface...%d-DIM: %d\n", (int)DIM, (int)N); 108 109 PetscCall(MatCreateFFT(PETSC_COMM_WORLD, DIM, dim, MATFFTW, &A)); 110 111 /* Create vectors that are compatible with parallel layout of A - must call MatCreateVecs()! */ 112 113 PetscCall(MatCreateVecsFFTW(A, &x, &y, &z)); 114 PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector")); 115 PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector")); 116 PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector")); 117 118 /* Set values of space vector x */ 119 PetscCall(VecSetRandom(x, rdm)); 120 121 if (view) PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 122 123 /* Apply FFTW_FORWARD and FFTW_BACKWARD */ 124 PetscCall(MatMult(A, x, y)); 125 if (view) PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD)); 126 127 PetscCall(MatMultTranspose(A, y, z)); 128 129 /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */ 130 a = 1.0 / (PetscReal)N; 131 PetscCall(VecScale(z, a)); 132 if (view) PetscCall(VecView(z, PETSC_VIEWER_STDOUT_WORLD)); 133 PetscCall(VecAXPY(z, -1.0, x)); 134 PetscCall(VecNorm(z, NORM_1, &enorm)); 135 if (enorm > 1.e-9 && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %e\n", enorm)); 136 137 PetscCall(VecDestroy(&x)); 138 PetscCall(VecDestroy(&y)); 139 PetscCall(VecDestroy(&z)); 140 PetscCall(MatDestroy(&A)); 141 142 PetscCall(PetscFree(dim)); 143 } 144 } 145 146 PetscCall(PetscRandomDestroy(&rdm)); 147 PetscCall(PetscFinalize()); 148 return 0; 149 } 150 151 /*TEST 152 153 build: 154 requires: !mpiuni fftw complex 155 156 test: 157 output_file: output/ex143.out 158 159 test: 160 suffix: 2 161 nsize: 3 162 output_file: output/ex143.out 163 164 TEST*/ 165