Actual source code: ex148.c
petsc-3.14.6 2021-03-30
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";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
9: PetscMPIInt rank,size;
10: PetscInt N0=50,N1=50,N=N0*N1;
11: PetscRandom rdm;
12: PetscReal enorm;
13: Vec x,y,z,input,output;
14: Mat A;
15: PetscInt DIM,dim[2];
16: PetscReal fac;
18: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
19: #if defined(PETSC_USE_COMPLEX)
20: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires real numbers");
21: #endif
23: MPI_Comm_size(PETSC_COMM_WORLD, &size);
24: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
26: /* Create and set PETSc vectors 'input' and 'output' */
27: PetscRandomCreate(PETSC_COMM_WORLD, &rdm);
28: PetscRandomSetFromOptions(rdm);
30: VecCreate(PETSC_COMM_WORLD,&input);
31: VecSetSizes(input,PETSC_DECIDE,N0*N1);
32: VecSetFromOptions(input);
33: VecSetRandom(input,rdm);
34: VecDuplicate(input,&output);
35: PetscObjectSetName((PetscObject)input, "Real space vector");
36: PetscObjectSetName((PetscObject)output, "Reconstructed vector");
38: /* Get FFTW vectors 'x', 'y' and 'z' */
39: DIM = 2;
40: dim[0] = N0; dim[1] = N1;
41: MatCreateFFT(PETSC_COMM_WORLD,DIM,dim,MATFFTW,&A);
42: MatCreateVecsFFTW(A,&x,&y,&z);
44: /* Scatter PETSc vector 'input' to FFTW vector 'x' */
45: VecScatterPetscToFFTW(A,input,x);
47: /* Apply forward FFT */
48: MatMult(A,x,y);
50: /* Apply backward FFT */
51: MatMultTranspose(A,y,z);
53: /* Scatter FFTW vector 'z' to PETSc vector 'output' */
54: VecScatterFFTWToPetsc(A,z,output);
56: /* Check accuracy */
57: fac = 1.0/(PetscReal)N;
58: VecScale(output,fac);
59: VecAXPY(output,-1.0,input);
60: VecNorm(output,NORM_1,&enorm);
61: if (enorm > 1.e-11 && !rank) {
62: PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %e\n",enorm);
63: }
65: /* Free spaces */
66: PetscRandomDestroy(&rdm);
67: VecDestroy(&input);
68: VecDestroy(&output);
69: VecDestroy(&x);
70: VecDestroy(&y);
71: VecDestroy(&z);
72: MatDestroy(&A);
74: PetscFinalize();
75: return ierr;
76: }
82: /*TEST
84: build:
85: requires: fftw !complex
87: test:
88: output_file: output/ex148.out
90: test:
91: suffix: 2
92: nsize: 3
93: output_file: output/ex148.out
95: TEST*/