Actual source code: ex228.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      1: static char help[] = "Test duplication/destruction of FFTW vecs \n\n";

      3: /*
      4:  Compiling the code:
      5:    This code uses the FFTW interface.
      6:    Use one of the options below to configure:
      7:    --with-fftw-dir=/.... or --download-fftw
      8:  Usage:
      9:    mpiexec -np <np> ./ex228
     10: */

     12:  #include <petscmat.h>
     13: int main(int argc,char **args)
     14: {
     15:   Mat            A;         /* FFT Matrix */
     16:   Vec            x,y,z;     /* Work vectors */
     17:   Vec            x1,y1,z1;  /* Duplicate vectors */
     18:   PetscInt       i,k;       /* for iterating over dimensions */
     19:   PetscRandom    rdm;       /* for creating random input */
     20:   PetscScalar    a;         /* used to scale output */
     21:   PetscReal      enorm;     /* norm for sanity check */
     23:   PetscInt       n=10,N=1;  /* FFT dimension params */
     24:   PetscInt       DIM,dim[5];/* FFT params */

     26:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
     27:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);

     29:   /* To create random input vector */
     30:   PetscRandomCreate(PETSC_COMM_SELF, &rdm);
     31:   PetscRandomSetFromOptions(rdm);

     33:   /* Iterate over dimensions, use PETSc-FFTW interface */
     34:   for (i=1; i<5; i++) {
     35:     DIM = i;
     36:     N = 1;
     37:     for (k=0; k<i; k++){dim[k] = n; N*=n;}

     39:     PetscPrintf(PETSC_COMM_WORLD, "\n %D dimensions: FFTW on vector of size %D \n",DIM,N);

     41:     /* create FFTW object */
     42:     MatCreateFFT(PETSC_COMM_SELF,DIM,dim,MATFFTW,&A);
     43:     /* create vectors of length N */
     44:     MatCreateVecsFFTW(A,&x,&y,&z);

     46:     PetscObjectSetName((PetscObject) x, "Real space vector");
     47:     PetscObjectSetName((PetscObject) y, "Frequency space vector");
     48:     PetscObjectSetName((PetscObject) z, "Reconstructed vector");

     50:     /* Test vector duplication*/
     51:     VecDuplicate(x,&x1);
     52:     VecDuplicate(y,&y1);
     53:     VecDuplicate(z,&z1);

     55:     /* Set values of space vector x, copy to duplicate */
     56:     VecSetRandom(x,rdm);
     57:     VecCopy(x,x1);

     59:     /* Apply FFTW_FORWARD and FFTW_BACKWARD */
     60:     MatMult(A,x,y);
     61:     MatMultTranspose(A,y,z);

     63:     /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */
     64:     MatMult(A,x1,y1);
     65:     MatMultTranspose(A,y1,z1);

     67:     /* Compare x and z1. FFTW computes an unnormalized DFT, thus z1 = N*x */
     68:     a    = 1.0/(PetscReal)N;
     69:     VecScale(z1,a);
     70:     VecAXPY(z1,-1.0,x);
     71:     VecNorm(z1,NORM_1,&enorm);
     72:     if (enorm > 1.e-9){PetscPrintf(PETSC_COMM_WORLD,"  Error norm of |x - z1| %g\n",enorm);}

     74:     /* free spaces */
     75:     VecDestroy(&x1);
     76:     VecDestroy(&y1);
     77:     VecDestroy(&z1);
     78:     VecDestroy(&x);
     79:     VecDestroy(&y);
     80:     VecDestroy(&z);
     81:     MatDestroy(&A);
     82:   }

     84:   PetscRandomDestroy(&rdm);
     85:   PetscFinalize();
     86:   return ierr;
     87: }

     89: /*TEST

     91:     build:
     92:       requires: fftw complex

     94:     test:
     95:       suffix: 2
     96:       nsize : 4
     97:       args: -mat_fftw_plannerflags FFTW_ESTIMATE -n 16

     99:     test:
    100:       suffix: 3
    101:       nsize : 2
    102:       args: -mat_fftw_plannerflags FFTW_MEASURE -n 12

    104:     test:
    105:       suffix: 4
    106:       nsize : 2
    107:       args: -mat_fftw_plannerflags FFTW_PATIENT -n 10

    109:     test:
    110:       suffix: 5
    111:       nsize : 1
    112:       args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE -n 5

    114: TEST*/