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