Actual source code: ex148.c
1: static char help[]="This program illustrates the use of PETSc-fftw interface for real 2D DFT\n";
2: #include <petscmat.h>
6: PetscInt main(PetscInt 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);
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");
37:
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: MatGetVecsFFTW(A,&x,&y,&z);
44: /* Scatter PETSc vector 'input' to FFTW vector 'x' */
45: InputTransformFFT(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: OutputTransformFFT(A,z,output);
56: /* Check accuracy */
57: fac = 1.0/(PetscReal)N;
58: VecScale(output,fac);
59: // VecView(input,PETSC_VIEWER_STDOUT_WORLD);
60: // VecView(output,PETSC_VIEWER_STDOUT_WORLD);
62: VecAXPY(output,-1.0,input);
63: VecNorm(output,NORM_1,&enorm);
64: if (enorm > 1.e-11 && !rank){
65: PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %e\n",enorm);
66: }
68: /* Free spaces */
69: PetscRandomDestroy(&rdm);
70: VecDestroy(&input);
71: VecDestroy(&output);
72: VecDestroy(&x);
73: VecDestroy(&y);
74: VecDestroy(&z);
75: MatDestroy(&A);
77: PetscFinalize();
78: return 0;
79: }