xref: /petsc/src/mat/tests/ex150.c (revision c4762a1b19cd2af06abeed90e8f9d34fb975dd94) 
1*c4762a1bSJed Brown static char help[]="This program illustrates the use of PETSc-fftw interface for real DFT\n";
2*c4762a1bSJed Brown #include <petscmat.h>
3*c4762a1bSJed Brown #include <fftw3-mpi.h>
4*c4762a1bSJed Brown 
5*c4762a1bSJed Brown extern PetscErrorCode InputTransformFFT(Mat,Vec,Vec);
6*c4762a1bSJed Brown extern PetscErrorCode OutputTransformFFT(Mat,Vec,Vec);
7*c4762a1bSJed Brown PetscInt main(PetscInt argc,char **args)
8*c4762a1bSJed Brown {
9*c4762a1bSJed Brown   PetscErrorCode ierr;
10*c4762a1bSJed Brown   PetscMPIInt    rank,size;
11*c4762a1bSJed Brown   PetscInt       N0=3,N1=3,N2=3,N3=3,N4=3,N=N0*N1*N2*N3;
12*c4762a1bSJed Brown   PetscRandom    rdm;
13*c4762a1bSJed Brown   PetscReal      enorm;
14*c4762a1bSJed Brown   Vec            x,y,z,input,output;
15*c4762a1bSJed Brown   Mat            A;
16*c4762a1bSJed Brown   PetscInt       DIM, dim[5],vsize;
17*c4762a1bSJed Brown   PetscReal      fac;
18*c4762a1bSJed Brown   PetscScalar    one=1;
19*c4762a1bSJed Brown 
20*c4762a1bSJed Brown   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
21*c4762a1bSJed Brown #if defined(PETSC_USE_COMPLEX)
22*c4762a1bSJed Brown   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires real numbers");
23*c4762a1bSJed Brown #endif
24*c4762a1bSJed Brown 
25*c4762a1bSJed Brown   ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);
26*c4762a1bSJed Brown   ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
27*c4762a1bSJed Brown 
28*c4762a1bSJed Brown   ierr = PetscRandomCreate(PETSC_COMM_WORLD, &rdm);CHKERRQ(ierr);
29*c4762a1bSJed Brown   ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr);
30*c4762a1bSJed Brown   ierr = VecCreate(PETSC_COMM_WORLD,&input);CHKERRQ(ierr);
31*c4762a1bSJed Brown   ierr = VecSetSizes(input,PETSC_DECIDE,N);CHKERRQ(ierr);
32*c4762a1bSJed Brown   ierr = VecSetFromOptions(input);CHKERRQ(ierr);
33*c4762a1bSJed Brown   ierr = VecSetRandom(input,rdm);CHKERRQ(ierr);
34*c4762a1bSJed Brown   ierr = VecAssemblyBegin(input);CHKERRQ(ierr);
35*c4762a1bSJed Brown   ierr = VecAssemblyEnd(input);CHKERRQ(ierr);
36*c4762a1bSJed Brown /*  ierr = VecView(input,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
37*c4762a1bSJed Brown   ierr = VecDuplicate(input,&output);
38*c4762a1bSJed Brown 
39*c4762a1bSJed Brown   DIM    = 4;
40*c4762a1bSJed Brown   dim[0] = N0; dim[1] = N1; dim[2] = N2; dim[3] = N3; dim[4] = N4;
41*c4762a1bSJed Brown 
42*c4762a1bSJed Brown   ierr = MatCreateFFT(PETSC_COMM_WORLD,DIM,dim,MATFFTW,&A);CHKERRQ(ierr);
43*c4762a1bSJed Brown   ierr = MatCreateVecs(A,&x,&y);CHKERRQ(ierr);
44*c4762a1bSJed Brown   ierr = MatCreateVecs(A,&z,NULL);CHKERRQ(ierr);
45*c4762a1bSJed Brown   ierr = VecGetSize(x,&vsize);CHKERRQ(ierr);
46*c4762a1bSJed Brown   printf("The vector size from the main routine is %d\n",vsize);
47*c4762a1bSJed Brown 
48*c4762a1bSJed Brown   ierr = InputTransformFFT(A,input,x);CHKERRQ(ierr);
49*c4762a1bSJed Brown 
50*c4762a1bSJed Brown   ierr = MatMult(A,x,y);CHKERRQ(ierr);
51*c4762a1bSJed Brown   ierr = VecAssemblyBegin(y);CHKERRQ(ierr);
52*c4762a1bSJed Brown   ierr = VecAssemblyEnd(y);CHKERRQ(ierr);
53*c4762a1bSJed Brown   ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
54*c4762a1bSJed Brown   ierr = MatMultTranspose(A,y,z);CHKERRQ(ierr);
55*c4762a1bSJed Brown   ierr = VecGetSize(z,&vsize);CHKERRQ(ierr);
56*c4762a1bSJed Brown   printf("The vector size of zfrom the main routine is %d\n",vsize);
57*c4762a1bSJed Brown 
58*c4762a1bSJed Brown   ierr = OutputTransformFFT(A,z,output);CHKERRQ(ierr);
59*c4762a1bSJed Brown 
60*c4762a1bSJed Brown   fac  = 1.0/(PetscReal)N;
61*c4762a1bSJed Brown   ierr = VecScale(output,fac);CHKERRQ(ierr);
62*c4762a1bSJed Brown 
63*c4762a1bSJed Brown   ierr = VecAssemblyBegin(input);CHKERRQ(ierr);
64*c4762a1bSJed Brown   ierr = VecAssemblyEnd(input);CHKERRQ(ierr);
65*c4762a1bSJed Brown   ierr = VecAssemblyBegin(output);CHKERRQ(ierr);
66*c4762a1bSJed Brown   ierr = VecAssemblyEnd(output);CHKERRQ(ierr);
67*c4762a1bSJed Brown 
68*c4762a1bSJed Brown   ierr = VecView(input,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
69*c4762a1bSJed Brown   ierr = VecView(output,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
70*c4762a1bSJed Brown 
71*c4762a1bSJed Brown   ierr = VecAXPY(output,-1.0,input);CHKERRQ(ierr);
72*c4762a1bSJed Brown   ierr = VecNorm(output,NORM_1,&enorm);CHKERRQ(ierr);
73*c4762a1bSJed Brown /*  if (enorm > 1.e-14) { */
74*c4762a1bSJed Brown   if (!rank) {
75*c4762a1bSJed Brown     ierr = PetscPrintf(PETSC_COMM_SELF,"  Error norm of |x - z| %e\n",enorm);CHKERRQ(ierr);
76*c4762a1bSJed Brown   }
77*c4762a1bSJed Brown /*      } */
78*c4762a1bSJed Brown 
79*c4762a1bSJed Brown 
80*c4762a1bSJed Brown 
81*c4762a1bSJed Brown 
82*c4762a1bSJed Brown /* ierr = MatCreateVecs(A,&z,NULL);CHKERRQ(ierr); */
83*c4762a1bSJed Brown /*  printf("Vector size from ex148 %d\n",vsize); */
84*c4762a1bSJed Brown /*  ierr = PetscObjectSetName((PetscObject) x, "Real space vector");CHKERRQ(ierr); */
85*c4762a1bSJed Brown /*      ierr = PetscObjectSetName((PetscObject) y, "Frequency space vector");CHKERRQ(ierr); */
86*c4762a1bSJed Brown /*      ierr = PetscObjectSetName((PetscObject) z, "Reconstructed vector");CHKERRQ(ierr); */
87*c4762a1bSJed Brown 
88*c4762a1bSJed Brown   ierr = VecDestroy(&output);CHKERRQ(ierr);
89*c4762a1bSJed Brown   ierr = VecDestroy(&input);CHKERRQ(ierr);
90*c4762a1bSJed Brown   ierr = VecDestroy(&x);CHKERRQ(ierr);
91*c4762a1bSJed Brown   ierr = VecDestroy(&y);CHKERRQ(ierr);
92*c4762a1bSJed Brown   ierr = VecDestroy(&z);CHKERRQ(ierr);
93*c4762a1bSJed Brown   ierr = MatDestroy(&A);CHKERRQ(ierr);
94*c4762a1bSJed Brown   ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr);
95*c4762a1bSJed Brown   ierr = PetscFinalize();
96*c4762a1bSJed Brown   return ierr;
97*c4762a1bSJed Brown }
98*c4762a1bSJed Brown 
99*c4762a1bSJed Brown 
100*c4762a1bSJed Brown 
101