1 static char help[] = "This program illustrates the use of PETSc-fftw interface for parallel real DFT\n";
2 #include <petscmat.h>
3 #include <fftw3-mpi.h>
4 /*extern PetscErrorCode MatCreateVecsFFT(Mat,Vec *,Vec *,Vec *);*/
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 PetscMPIInt rank, size;
8 PetscInt N0 = 4096, N1 = 4096, N2 = 256, N3 = 10, N4 = 10, N = N0 * N1;
9 PetscRandom rdm;
10 PetscReal enorm;
11 Vec x, y, z, input, output;
12 Mat A;
13 PetscInt DIM, dim[5], vsize, row, col;
14 PetscReal fac;
15
16 PetscFunctionBeginUser;
17 PetscCall(PetscInitialize(&argc, &args, NULL, help));
18 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
19 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
20
21 #if defined(PETSC_USE_COMPLEX)
22 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Example for Real DFT. Your current data type is complex!");
23 #endif
24 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
25 PetscCall(PetscRandomSetFromOptions(rdm));
26 PetscCall(VecCreate(PETSC_COMM_WORLD, &input));
27 PetscCall(VecSetSizes(input, PETSC_DECIDE, N));
28 PetscCall(VecSetFromOptions(input));
29 PetscCall(VecSetRandom(input, rdm));
30 PetscCall(VecDuplicate(input, &output));
31
32 DIM = 2;
33 dim[0] = N0;
34 dim[1] = N1;
35 dim[2] = N2;
36 dim[3] = N3;
37 dim[4] = N4;
38 PetscCall(MatCreateFFT(PETSC_COMM_WORLD, DIM, dim, MATFFTW, &A));
39 PetscCall(MatGetLocalSize(A, &row, &col));
40 printf("The Matrix size is %d and %d from process %d\n", row, col, rank);
41 PetscCall(MatCreateVecsFFTW(A, &x, &y, &z));
42
43 PetscCall(VecGetSize(x, &vsize));
44
45 PetscCall(VecGetSize(z, &vsize));
46 printf("The vector size of output from the main routine is %d\n", vsize);
47
48 PetscCall(VecScatterPetscToFFTW(A, input, x));
49 /*PetscCall(VecDestroy(&input));*/
50 PetscCall(MatMult(A, x, y));
51 PetscCall(MatMultTranspose(A, y, z));
52 PetscCall(VecScatterFFTWToPetsc(A, z, output));
53 /*PetscCall(VecDestroy(&z));*/
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) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %e\n", enorm));
68
69 PetscCall(VecDestroy(&x));
70 PetscCall(VecDestroy(&y));
71 PetscCall(VecDestroy(&z));
72 PetscCall(VecDestroy(&output));
73 PetscCall(VecDestroy(&input));
74 PetscCall(MatDestroy(&A));
75 PetscCall(PetscRandomDestroy(&rdm));
76 PetscCall(PetscFinalize());
77 return 0;
78 }
79