1: static char help[] = "Test sequential FFTW interface \n\n";

  3: /*
  4:   Compiling the code:
  5:       This code uses the complex numbers version of PETSc, so configure
  6:       must be run to enable this

  8: */

 10:  #include petscmat.h
 13: PetscInt main(PetscInt argc,char **args)
 14: {
 15:   Mat            A;
 17:   PetscMPIInt    size;
 18:   PetscInt       n = 10,N,ndim=4,dim[4],DIM,i;
 19:   Vec            x,y,z;
 20:   PetscScalar    s;
 21:   PetscRandom    rdm;
 22:   PetscReal      enorm;

 24:   PetscInitialize(&argc,&args,(char *)0,help);
 25: #if !defined(PETSC_USE_COMPLEX)
 26:   SETERRQ(1,"This example requires complex numbers");
 27: #endif
 28:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 29:   if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");

 31:   for (DIM=0; DIM<ndim; DIM++){
 32:     dim[DIM] = n;  /* size of transformation in DIM-dimension */
 33:   }
 34:   PetscRandomCreate(PETSC_COMM_SELF,&rdm);
 35:   PetscRandomSetFromOptions(rdm);

 37:   for (DIM=1; DIM<5; DIM++){
 38:     /* create vectors of length N=n^DIM */
 39:     N = 1; for (i=0; i<DIM; i++) N *= dim[i];
 40:     PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n",DIM,N);
 41:     VecCreateSeq(PETSC_COMM_SELF,N,&x);
 42:     VecDuplicate(x,&y);
 43:     VecDuplicate(x,&z);
 44:     VecSetRandom(x,rdm);

 46:     /* create FFTW object */
 47:     MatCreateSeqFFTW(PETSC_COMM_SELF,DIM,dim,&A);

 49:     /* apply FFTW_FORWARD several times, so the fftw_plan can be reused on different vectors */
 50:     MatMult(A,x,z);
 51:     for (i=0; i<3; i++){
 52:       MatMult(A,x,y);

 54:       /* apply FFTW_BACKWARD several times */
 55:       MatMultTranspose(A,y,z);
 56:     }
 57: 
 58:     /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
 59:     s = 1.0/(PetscReal)N;
 60:     VecScale(z,s);
 61:     VecAXPY(z,-1.0,x);
 62:     VecNorm(z,NORM_1,&enorm);
 63:     PetscPrintf(PETSC_COMM_SELF,"  Error norm of |x - z| = %A\n",enorm);

 65:     /* free spaces */
 66:     VecDestroy(x);
 67:     VecDestroy(y);
 68:     VecDestroy(z);
 69:     MatDestroy(A);
 70:   }
 71:   PetscRandomDestroy(rdm);
 72:   PetscFinalize();
 73:   return 0;
 74: }