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