1 static char help[] = "This program illustrates the use of PETSc-fftw interface for real DFT\n";
2 #include <petscmat.h>
3 #include <fftw3-mpi.h>
4
5 extern PetscErrorCode InputTransformFFT(Mat, Vec, Vec);
6 extern PetscErrorCode OutputTransformFFT(Mat, Vec, Vec);
main(int argc,char ** args)7 int main(int argc, char **args)
8 {
9 PetscMPIInt rank, size;
10 PetscInt N0 = 3, N1 = 3, N2 = 3, N = N0 * N1 * N2;
11 PetscRandom rdm;
12 PetscScalar a;
13 PetscReal enorm;
14 Vec x, y, z, input, output;
15 PetscBool view = PETSC_FALSE, use_interface = PETSC_TRUE;
16 Mat A;
17 PetscInt DIM, dim[3], vsize;
18 PetscReal fac;
19
20 PetscFunctionBeginUser;
21 PetscCall(PetscInitialize(&argc, &args, NULL, help));
22 PetscCheck(!PetscDefined(USE_COMPLEX), PETSC_COMM_WORLD, PETSC_ERR_SUP, "This example requires real numbers");
23
24 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
25 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
26
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, N));
32 PetscCall(VecSetFromOptions(input));
33 PetscCall(VecSetRandom(input, rdm));
34 PetscCall(VecDuplicate(input, &output));
35 /* PetscCall(VecGetSize(input,&vsize)); */
36 /* printf("Size of the input Vector is %d\n",vsize); */
37
38 DIM = 3;
39 dim[0] = N0;
40 dim[1] = N1;
41 dim[2] = N2;
42
43 PetscCall(MatCreateFFT(PETSC_COMM_WORLD, DIM, dim, MATFFTW, &A));
44 PetscCall(MatCreateVecs(A, &x, &y));
45 PetscCall(MatCreateVecs(A, &z, NULL));
46 PetscCall(VecGetSize(y, &vsize));
47 printf("The vector size from the main routine is %d\n", vsize);
48
49 PetscCall(InputTransformFFT(A, input, x));
50 PetscCall(MatMult(A, x, y));
51 PetscCall(MatMultTranspose(A, y, z));
52 PetscCall(OutputTransformFFT(A, z, output));
53
54 fac = 1.0 / (PetscReal)N;
55 PetscCall(VecScale(output, fac));
56
57 PetscCall(VecAssemblyBegin(input));
58 PetscCall(VecAssemblyEnd(input));
59 PetscCall(VecAssemblyBegin(output));
60 PetscCall(VecAssemblyEnd(output));
61
62 PetscCall(VecView(input, PETSC_VIEWER_STDOUT_WORLD));
63 PetscCall(VecView(output, PETSC_VIEWER_STDOUT_WORLD));
64
65 PetscCall(VecAXPY(output, -1.0, input));
66 PetscCall(VecNorm(output, NORM_1, &enorm));
67 /* if (enorm > 1.e-14) { */
68 if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %e\n", enorm));
69 /* } */
70
71 /* PetscCall(MatCreateVecs(A,&z,NULL)); */
72 /* printf("Vector size from ex148 %d\n",vsize); */
73 /* PetscCall(PetscObjectSetName((PetscObject) x, "Real space vector")); */
74 /* PetscCall(PetscObjectSetName((PetscObject) y, "Frequency space vector")); */
75 /* PetscCall(PetscObjectSetName((PetscObject) z, "Reconstructed vector")); */
76
77 PetscCall(PetscFinalize());
78 return 0;
79 }
80