1*c4762a1bSJed Brown static char help[] = "Test duplication/destruction of FFTW vecs \n\n"; 2*c4762a1bSJed Brown 3*c4762a1bSJed Brown /* 4*c4762a1bSJed Brown Compiling the code: 5*c4762a1bSJed Brown This code uses the FFTW interface. 6*c4762a1bSJed Brown Use one of the options below to configure: 7*c4762a1bSJed Brown --with-fftw-dir=/.... or --download-fftw 8*c4762a1bSJed Brown Usage: 9*c4762a1bSJed Brown mpiexec -np <np> ./ex228 10*c4762a1bSJed Brown */ 11*c4762a1bSJed Brown 12*c4762a1bSJed Brown #include <petscmat.h> 13*c4762a1bSJed Brown int main(int argc,char **args) 14*c4762a1bSJed Brown { 15*c4762a1bSJed Brown Mat A; /* FFT Matrix */ 16*c4762a1bSJed Brown Vec x,y,z; /* Work vectors */ 17*c4762a1bSJed Brown Vec x1,y1,z1; /* Duplicate vectors */ 18*c4762a1bSJed Brown PetscInt i,k; /* for iterating over dimensions */ 19*c4762a1bSJed Brown PetscRandom rdm; /* for creating random input */ 20*c4762a1bSJed Brown PetscScalar a; /* used to scale output */ 21*c4762a1bSJed Brown PetscReal enorm; /* norm for sanity check */ 22*c4762a1bSJed Brown PetscErrorCode ierr; /* to catch bugs, if any */ 23*c4762a1bSJed Brown PetscInt n=10,N=1; /* FFT dimension params */ 24*c4762a1bSJed Brown PetscInt DIM,dim[5];/* FFT params */ 25*c4762a1bSJed Brown 26*c4762a1bSJed Brown ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 27*c4762a1bSJed Brown ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 28*c4762a1bSJed Brown 29*c4762a1bSJed Brown /* To create random input vector */ 30*c4762a1bSJed Brown ierr = PetscRandomCreate(PETSC_COMM_SELF, &rdm);CHKERRQ(ierr); 31*c4762a1bSJed Brown ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr); 32*c4762a1bSJed Brown 33*c4762a1bSJed Brown /* Iterate over dimensions, use PETSc-FFTW interface */ 34*c4762a1bSJed Brown for (i=1; i<5; i++) { 35*c4762a1bSJed Brown DIM = i; 36*c4762a1bSJed Brown N = 1; 37*c4762a1bSJed Brown for (k=0; k<i; k++){dim[k] = n; N*=n;} 38*c4762a1bSJed Brown 39*c4762a1bSJed Brown ierr = PetscPrintf(PETSC_COMM_WORLD, "\n %D dimensions: FFTW on vector of size %D \n",DIM,N);CHKERRQ(ierr); 40*c4762a1bSJed Brown 41*c4762a1bSJed Brown /* create FFTW object */ 42*c4762a1bSJed Brown ierr = MatCreateFFT(PETSC_COMM_SELF,DIM,dim,MATFFTW,&A);CHKERRQ(ierr); 43*c4762a1bSJed Brown /* create vectors of length N */ 44*c4762a1bSJed Brown ierr = MatCreateVecsFFTW(A,&x,&y,&z);CHKERRQ(ierr); 45*c4762a1bSJed Brown 46*c4762a1bSJed Brown ierr = PetscObjectSetName((PetscObject) x, "Real space vector");CHKERRQ(ierr); 47*c4762a1bSJed Brown ierr = PetscObjectSetName((PetscObject) y, "Frequency space vector");CHKERRQ(ierr); 48*c4762a1bSJed Brown ierr = PetscObjectSetName((PetscObject) z, "Reconstructed vector");CHKERRQ(ierr); 49*c4762a1bSJed Brown 50*c4762a1bSJed Brown /* Test vector duplication*/ 51*c4762a1bSJed Brown ierr = VecDuplicate(x,&x1);CHKERRQ(ierr); 52*c4762a1bSJed Brown ierr = VecDuplicate(y,&y1);CHKERRQ(ierr); 53*c4762a1bSJed Brown ierr = VecDuplicate(z,&z1);CHKERRQ(ierr); 54*c4762a1bSJed Brown 55*c4762a1bSJed Brown /* Set values of space vector x, copy to duplicate */ 56*c4762a1bSJed Brown ierr = VecSetRandom(x,rdm);CHKERRQ(ierr); 57*c4762a1bSJed Brown ierr = VecCopy(x,x1);CHKERRQ(ierr); 58*c4762a1bSJed Brown 59*c4762a1bSJed Brown /* Apply FFTW_FORWARD and FFTW_BACKWARD */ 60*c4762a1bSJed Brown ierr = MatMult(A,x,y);CHKERRQ(ierr); 61*c4762a1bSJed Brown ierr = MatMultTranspose(A,y,z);CHKERRQ(ierr); 62*c4762a1bSJed Brown 63*c4762a1bSJed Brown /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */ 64*c4762a1bSJed Brown ierr = MatMult(A,x1,y1);CHKERRQ(ierr); 65*c4762a1bSJed Brown ierr = MatMultTranspose(A,y1,z1);CHKERRQ(ierr); 66*c4762a1bSJed Brown 67*c4762a1bSJed Brown /* Compare x and z1. FFTW computes an unnormalized DFT, thus z1 = N*x */ 68*c4762a1bSJed Brown a = 1.0/(PetscReal)N; 69*c4762a1bSJed Brown ierr = VecScale(z1,a);CHKERRQ(ierr); 70*c4762a1bSJed Brown ierr = VecAXPY(z1,-1.0,x);CHKERRQ(ierr); 71*c4762a1bSJed Brown ierr = VecNorm(z1,NORM_1,&enorm);CHKERRQ(ierr); 72*c4762a1bSJed Brown if (enorm > 1.e-9){ierr = PetscPrintf(PETSC_COMM_WORLD," Error norm of |x - z1| %g\n",enorm);CHKERRQ(ierr);} 73*c4762a1bSJed Brown 74*c4762a1bSJed Brown /* free spaces */ 75*c4762a1bSJed Brown ierr = VecDestroy(&x1);CHKERRQ(ierr); 76*c4762a1bSJed Brown ierr = VecDestroy(&y1);CHKERRQ(ierr); 77*c4762a1bSJed Brown ierr = VecDestroy(&z1);CHKERRQ(ierr); 78*c4762a1bSJed Brown ierr = VecDestroy(&x);CHKERRQ(ierr); 79*c4762a1bSJed Brown ierr = VecDestroy(&y);CHKERRQ(ierr); 80*c4762a1bSJed Brown ierr = VecDestroy(&z);CHKERRQ(ierr); 81*c4762a1bSJed Brown ierr = MatDestroy(&A);CHKERRQ(ierr); 82*c4762a1bSJed Brown } 83*c4762a1bSJed Brown 84*c4762a1bSJed Brown ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr); 85*c4762a1bSJed Brown ierr = PetscFinalize(); 86*c4762a1bSJed Brown return ierr; 87*c4762a1bSJed Brown } 88*c4762a1bSJed Brown 89*c4762a1bSJed Brown /*TEST 90*c4762a1bSJed Brown 91*c4762a1bSJed Brown build: 92*c4762a1bSJed Brown requires: fftw complex 93*c4762a1bSJed Brown 94*c4762a1bSJed Brown test: 95*c4762a1bSJed Brown suffix: 2 96*c4762a1bSJed Brown nsize : 4 97*c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_ESTIMATE -n 16 98*c4762a1bSJed Brown 99*c4762a1bSJed Brown test: 100*c4762a1bSJed Brown suffix: 3 101*c4762a1bSJed Brown nsize : 2 102*c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_MEASURE -n 12 103*c4762a1bSJed Brown 104*c4762a1bSJed Brown test: 105*c4762a1bSJed Brown suffix: 4 106*c4762a1bSJed Brown nsize : 2 107*c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_PATIENT -n 10 108*c4762a1bSJed Brown 109*c4762a1bSJed Brown test: 110*c4762a1bSJed Brown suffix: 5 111*c4762a1bSJed Brown nsize : 1 112*c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE -n 5 113*c4762a1bSJed Brown 114*c4762a1bSJed Brown TEST*/ 115