1 static char help[] = "Test duplication/destruction of FFTW vecs \n\n"; 2 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 */ 11 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 */ 22 PetscErrorCode ierr; /* to catch bugs, if any */ 23 PetscInt n=10,N=1; /* FFT dimension params */ 24 PetscInt DIM,dim[5];/* FFT params */ 25 26 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 27 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 28 29 /* To create random input vector */ 30 ierr = PetscRandomCreate(PETSC_COMM_SELF, &rdm);CHKERRQ(ierr); 31 ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr); 32 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;} 38 39 ierr = PetscPrintf(PETSC_COMM_WORLD, "\n %D dimensions: FFTW on vector of size %D \n",DIM,N);CHKERRQ(ierr); 40 41 /* create FFTW object */ 42 ierr = MatCreateFFT(PETSC_COMM_SELF,DIM,dim,MATFFTW,&A);CHKERRQ(ierr); 43 /* create vectors of length N */ 44 ierr = MatCreateVecsFFTW(A,&x,&y,&z);CHKERRQ(ierr); 45 46 ierr = PetscObjectSetName((PetscObject) x, "Real space vector");CHKERRQ(ierr); 47 ierr = PetscObjectSetName((PetscObject) y, "Frequency space vector");CHKERRQ(ierr); 48 ierr = PetscObjectSetName((PetscObject) z, "Reconstructed vector");CHKERRQ(ierr); 49 50 /* Test vector duplication*/ 51 ierr = VecDuplicate(x,&x1);CHKERRQ(ierr); 52 ierr = VecDuplicate(y,&y1);CHKERRQ(ierr); 53 ierr = VecDuplicate(z,&z1);CHKERRQ(ierr); 54 55 /* Set values of space vector x, copy to duplicate */ 56 ierr = VecSetRandom(x,rdm);CHKERRQ(ierr); 57 ierr = VecCopy(x,x1);CHKERRQ(ierr); 58 59 /* Apply FFTW_FORWARD and FFTW_BACKWARD */ 60 ierr = MatMult(A,x,y);CHKERRQ(ierr); 61 ierr = MatMultTranspose(A,y,z);CHKERRQ(ierr); 62 63 /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */ 64 ierr = MatMult(A,x1,y1);CHKERRQ(ierr); 65 ierr = MatMultTranspose(A,y1,z1);CHKERRQ(ierr); 66 67 /* Compare x and z1. FFTW computes an unnormalized DFT, thus z1 = N*x */ 68 a = 1.0/(PetscReal)N; 69 ierr = VecScale(z1,a);CHKERRQ(ierr); 70 ierr = VecAXPY(z1,-1.0,x);CHKERRQ(ierr); 71 ierr = VecNorm(z1,NORM_1,&enorm);CHKERRQ(ierr); 72 if (enorm > 1.e-9){ierr = PetscPrintf(PETSC_COMM_WORLD," Error norm of |x - z1| %g\n",enorm);CHKERRQ(ierr);} 73 74 /* free spaces */ 75 ierr = VecDestroy(&x1);CHKERRQ(ierr); 76 ierr = VecDestroy(&y1);CHKERRQ(ierr); 77 ierr = VecDestroy(&z1);CHKERRQ(ierr); 78 ierr = VecDestroy(&x);CHKERRQ(ierr); 79 ierr = VecDestroy(&y);CHKERRQ(ierr); 80 ierr = VecDestroy(&z);CHKERRQ(ierr); 81 ierr = MatDestroy(&A);CHKERRQ(ierr); 82 } 83 84 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr); 85 ierr = PetscFinalize(); 86 return ierr; 87 } 88 89 /*TEST 90 91 build: 92 requires: fftw complex 93 94 test: 95 suffix: 2 96 nsize : 4 97 args: -mat_fftw_plannerflags FFTW_ESTIMATE -n 16 98 99 test: 100 suffix: 3 101 nsize : 2 102 args: -mat_fftw_plannerflags FFTW_MEASURE -n 12 103 104 test: 105 suffix: 4 106 nsize : 2 107 args: -mat_fftw_plannerflags FFTW_PATIENT -n 10 108 109 test: 110 suffix: 5 111 nsize : 1 112 args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE -n 5 113 114 TEST*/ 115