1 static char help[] = "This program illustrates the use of PETSc-fftw interface for real 2D DFT.\n\
2 See ~petsc/src/mat/tests/ex158.c for general cases. \n\n";
3
4 #include <petscmat.h>
5
main(int argc,char ** args)6 int main(int argc, char **args)
7 {
8 PetscMPIInt rank, size;
9 PetscInt N0 = 50, N1 = 50, N = N0 * N1;
10 PetscRandom rdm;
11 PetscReal enorm;
12 Vec x, y, z, input, output;
13 Mat A;
14 PetscInt DIM, dim[2];
15 PetscReal fac;
16
17 PetscFunctionBeginUser;
18 PetscCall(PetscInitialize(&argc, &args, NULL, help));
19 #if defined(PETSC_USE_COMPLEX)
20 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This example requires real numbers");
21 #endif
22
23 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
24 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
25
26 /* Create and set PETSc vectors 'input' and 'output' */
27 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
28 PetscCall(PetscRandomSetFromOptions(rdm));
29
30 PetscCall(VecCreate(PETSC_COMM_WORLD, &input));
31 PetscCall(VecSetSizes(input, PETSC_DECIDE, N0 * N1));
32 PetscCall(VecSetFromOptions(input));
33 PetscCall(VecSetRandom(input, rdm));
34 PetscCall(VecDuplicate(input, &output));
35 PetscCall(PetscObjectSetName((PetscObject)input, "Real space vector"));
36 PetscCall(PetscObjectSetName((PetscObject)output, "Reconstructed vector"));
37
38 /* Get FFTW vectors 'x', 'y' and 'z' */
39 DIM = 2;
40 dim[0] = N0;
41 dim[1] = N1;
42 PetscCall(MatCreateFFT(PETSC_COMM_WORLD, DIM, dim, MATFFTW, &A));
43 PetscCall(MatCreateVecsFFTW(A, &x, &y, &z));
44
45 /* Scatter PETSc vector 'input' to FFTW vector 'x' */
46 PetscCall(VecScatterPetscToFFTW(A, input, x));
47
48 /* Apply forward FFT */
49 PetscCall(MatMult(A, x, y));
50
51 /* Apply backward FFT */
52 PetscCall(MatMultTranspose(A, y, z));
53
54 /* Scatter FFTW vector 'z' to PETSc vector 'output' */
55 PetscCall(VecScatterFFTWToPetsc(A, z, output));
56
57 /* Check accuracy */
58 fac = 1.0 / (PetscReal)N;
59 PetscCall(VecScale(output, fac));
60 PetscCall(VecAXPY(output, -1.0, input));
61 PetscCall(VecNorm(output, NORM_1, &enorm));
62 if (enorm > 1.e-11 && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %e\n", enorm));
63
64 /* Free spaces */
65 PetscCall(PetscRandomDestroy(&rdm));
66 PetscCall(VecDestroy(&input));
67 PetscCall(VecDestroy(&output));
68 PetscCall(VecDestroy(&x));
69 PetscCall(VecDestroy(&y));
70 PetscCall(VecDestroy(&z));
71 PetscCall(MatDestroy(&A));
72
73 PetscCall(PetscFinalize());
74 return 0;
75 }
76
77 /*TEST
78
79 build:
80 requires: fftw !complex
81
82 test:
83 output_file: output/empty.out
84
85 test:
86 suffix: 2
87 nsize: 3
88 output_file: output/empty.out
89
90 TEST*/
91