1 static char help[] = "Test sequential FFTW interface \n\n"; 2 3 /* 4 Compiling the code: 5 This code uses the complex numbers version of PETSc, so configure 6 must be run to enable this 7 8 */ 9 10 #include <petscmat.h> 11 int main(int argc,char **args) 12 { 13 typedef enum {RANDOM, CONSTANT, TANH, NUM_FUNCS} FuncType; 14 const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"}; 15 Mat A; 16 PetscMPIInt size; 17 PetscInt n = 10,N,ndim=4,dim[4],DIM,i; 18 Vec x,y,z; 19 PetscScalar s; 20 PetscRandom rdm; 21 PetscReal enorm, tol = PETSC_SMALL; 22 PetscInt func; 23 FuncType function = RANDOM; 24 PetscBool view = PETSC_FALSE; 25 PetscErrorCode ierr; 26 27 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 28 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRMPI(ierr); 29 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This is a uniprocessor example only!"); 30 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex112");CHKERRQ(ierr); 31 ierr = PetscOptionsEList("-function", "Function type", "ex112", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);CHKERRQ(ierr); 32 ierr = PetscOptionsBool("-vec_view_draw", "View the functions", "ex112", view, &view, NULL);CHKERRQ(ierr); 33 function = (FuncType) func; 34 ierr = PetscOptionsEnd();CHKERRQ(ierr); 35 36 for (DIM = 0; DIM < ndim; DIM++) { 37 dim[DIM] = n; /* size of transformation in DIM-dimension */ 38 } 39 ierr = PetscRandomCreate(PETSC_COMM_SELF, &rdm);CHKERRQ(ierr); 40 ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr); 41 42 for (DIM = 1; DIM < 5; DIM++) { 43 for (i = 0, N = 1; i < DIM; i++) N *= dim[i]; 44 ierr = PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n",DIM,N);CHKERRQ(ierr); 45 46 /* create FFTW object */ 47 ierr = MatCreateFFT(PETSC_COMM_SELF,DIM,dim,MATFFTW,&A);CHKERRQ(ierr); 48 49 /* create vectors of length N=n^DIM */ 50 ierr = MatCreateVecs(A,&x,&y);CHKERRQ(ierr); 51 ierr = MatCreateVecs(A,&z,NULL);CHKERRQ(ierr); 52 ierr = PetscObjectSetName((PetscObject) x, "Real space vector");CHKERRQ(ierr); 53 ierr = PetscObjectSetName((PetscObject) y, "Frequency space vector");CHKERRQ(ierr); 54 ierr = PetscObjectSetName((PetscObject) z, "Reconstructed vector");CHKERRQ(ierr); 55 56 /* set values of space vector x */ 57 if (function == RANDOM) { 58 ierr = VecSetRandom(x, rdm);CHKERRQ(ierr); 59 } else if (function == CONSTANT) { 60 ierr = VecSet(x, 1.0);CHKERRQ(ierr); 61 } else if (function == TANH) { 62 PetscScalar *a; 63 ierr = VecGetArray(x, &a);CHKERRQ(ierr); 64 for (i = 0; i < N; ++i) { 65 a[i] = tanh((i - N/2.0)*(10.0/N)); 66 } 67 ierr = VecRestoreArray(x, &a);CHKERRQ(ierr); 68 } 69 if (view) {ierr = VecView(x, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);} 70 71 /* apply FFTW_FORWARD and FFTW_BACKWARD several times on same x, y, and z */ 72 for (i=0; i<3; i++) { 73 ierr = MatMult(A,x,y);CHKERRQ(ierr); 74 if (view && i == 0) {ierr = VecView(y, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);} 75 ierr = MatMultTranspose(A,y,z);CHKERRQ(ierr); 76 77 /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */ 78 s = 1.0/(PetscReal)N; 79 ierr = VecScale(z,s);CHKERRQ(ierr); 80 if (view && i == 0) {ierr = VecView(z, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);} 81 ierr = VecAXPY(z,-1.0,x);CHKERRQ(ierr); 82 ierr = VecNorm(z,NORM_1,&enorm);CHKERRQ(ierr); 83 if (enorm > tol) { 84 ierr = PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %g\n",(double)enorm);CHKERRQ(ierr); 85 } 86 } 87 88 /* apply FFTW_FORWARD and FFTW_BACKWARD several times on different x */ 89 for (i=0; i<3; i++) { 90 ierr = VecDestroy(&x);CHKERRQ(ierr); 91 ierr = VecCreateSeq(PETSC_COMM_SELF,N,&x);CHKERRQ(ierr); 92 ierr = VecSetRandom(x, rdm);CHKERRQ(ierr); 93 94 ierr = MatMult(A,x,y);CHKERRQ(ierr); 95 ierr = MatMultTranspose(A,y,z);CHKERRQ(ierr); 96 97 /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */ 98 s = 1.0/(PetscReal)N; 99 ierr = VecScale(z,s);CHKERRQ(ierr); 100 if (view && i == 0) {ierr = VecView(z, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);} 101 ierr = VecAXPY(z,-1.0,x);CHKERRQ(ierr); 102 ierr = VecNorm(z,NORM_1,&enorm);CHKERRQ(ierr); 103 if (enorm > tol) { 104 ierr = PetscPrintf(PETSC_COMM_SELF," Error norm of new |x - z| %g\n",(double)enorm);CHKERRQ(ierr); 105 } 106 } 107 108 /* free spaces */ 109 ierr = VecDestroy(&x);CHKERRQ(ierr); 110 ierr = VecDestroy(&y);CHKERRQ(ierr); 111 ierr = VecDestroy(&z);CHKERRQ(ierr); 112 ierr = MatDestroy(&A);CHKERRQ(ierr); 113 } 114 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr); 115 ierr = PetscFinalize(); 116 return ierr; 117 } 118 119 120 /*TEST 121 122 build: 123 requires: fftw complex 124 125 test: 126 args: -mat_fftw_plannerflags FFTW_ESTIMATE 127 output_file: output/ex112.out 128 129 test: 130 suffix: 2 131 args: -mat_fftw_plannerflags FFTW_MEASURE 132 output_file: output/ex112.out 133 requires: !define(PETSC_USE_CXXCOMPLEX) 134 135 test: 136 suffix: 3 137 args: -mat_fftw_plannerflags FFTW_PATIENT 138 output_file: output/ex112.out 139 140 test: 141 suffix: 4 142 args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE 143 output_file: output/ex112.out 144 145 TEST*/ 146