Actual source code: ex157.c
petsc-3.3-p6 2013-02-11
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>
6: PetscInt main(PetscInt argc,char **args)
7: {
8: PetscErrorCode ierr;
9: PetscMPIInt rank,size;
10: PetscInt N0=2048,N1=2048,N2=3,N3=5,N4=5,N=N0*N1;
11: PetscRandom rdm;
12: PetscReal enorm;
13: Vec x,y,z,input,output;
14: Mat A;
15: PetscInt DIM, dim[5],vsize;
16: PetscReal fac;
17: PetscScalar one=1,two=2,three=3;
19: PetscInitialize(&argc,&args,(char *)0,help);
20: MPI_Comm_size(PETSC_COMM_WORLD, &size);
21: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
22:
23: //#if !defined(PETSC_USE_COMPLEX)
24: // SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "Example for Complex DFT. Your current data type is real");
25: //#endif
26: PetscRandomCreate(PETSC_COMM_WORLD, &rdm);
27: PetscRandomSetFromOptions(rdm);
28: VecCreate(PETSC_COMM_WORLD,&input);
29: VecSetSizes(input,PETSC_DECIDE,N);
30: VecSetFromOptions(input);
31: // VecSet(input,one);
32: // VecSetValue(input,1,two,INSERT_VALUES);
33: // VecSetValue(input,2,three,INSERT_VALUES);
34: // VecSetValue(input,3,three,INSERT_VALUES);
35: VecSetRandom(input,rdm);
36: // VecSetRandom(input,rdm);
37: // VecSetRandom(input,rdm);
38: VecDuplicate(input,&output);
39:
40: DIM = 2; dim[0] = N0; dim[1] = N1; dim[2] = N2; dim[3] = N3; dim[4] = N4;
41: MatCreateFFT(PETSC_COMM_WORLD,DIM,dim,MATFFTW,&A);
42: MatGetVecsFFTW(A,&x,&y,&z);
43: // MatGetVecs(A,&x,&y);
44: // MatGetVecs(A,&z,PETSC_NULL);
46: VecGetSize(x,&vsize);
47: printf("The vector size of input from the main routine is %d\n",vsize);
49: VecGetSize(z,&vsize);
50: printf("The vector size of output from the main routine is %d\n",vsize);
52: InputTransformFFT(A,input,x);
54: MatMult(A,x,y);
55: VecAssemblyBegin(y);
56: VecAssemblyEnd(y);
57: VecView(y,PETSC_VIEWER_STDOUT_WORLD);
58:
59: MatMultTranspose(A,y,z);
61: OutputTransformFFT(A,z,output);
62: fac = 1.0/(PetscReal)N;
63: VecScale(output,fac);
64:
65: VecAssemblyBegin(input);
66: VecAssemblyEnd(input);
67: VecAssemblyBegin(output);
68: VecAssemblyEnd(output);
70: // VecView(input,PETSC_VIEWER_STDOUT_WORLD);
71: // VecView(output,PETSC_VIEWER_STDOUT_WORLD);
73: VecAXPY(output,-1.0,input);
74: VecNorm(output,NORM_1,&enorm);
75: // if (enorm > 1.e-14){
76: PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %e\n",enorm);
77: // }
79: VecDestroy(&output);
80: VecDestroy(&input);
81: VecDestroy(&x);
82: VecDestroy(&y);
83: VecDestroy(&z);
84: MatDestroy(&A);
85: PetscRandomDestroy(&rdm);
86: PetscFinalize();
87: return 0;
89: }