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