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