1 static char help[] = "Test sequential FFTW convolution\n\n"; 2 3 /* 4 Compiling the code: 5 This code uses the complex numbers, so configure must be given --with-scalar-type=complex to enable this 6 */ 7 8 #include <petscmat.h> 9 10 int main(int argc,char **args) 11 { 12 typedef enum {RANDOM, CONSTANT, TANH, NUM_FUNCS} FuncType; 13 const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"}; 14 Mat A; 15 PetscMPIInt size; 16 PetscInt n = 10,N,ndim=4,dim[4],DIM,i,j; 17 Vec w,x,y1,y2,z1,z2; 18 PetscScalar *a, *a2, *a3; 19 PetscScalar s; 20 PetscRandom rdm; 21 PetscReal enorm; 22 PetscInt func = 0; 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", "ex121", 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 /* create vectors of length N=n^DIM */ 43 for (i = 0, N = 1; i < DIM; i++) N *= dim[i]; 44 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n",DIM,N)); 45 PetscCall(VecCreateSeq(PETSC_COMM_SELF,N,&x)); 46 PetscCall(PetscObjectSetName((PetscObject) x, "Real space vector")); 47 PetscCall(VecDuplicate(x,&w)); 48 PetscCall(PetscObjectSetName((PetscObject) w, "Window vector")); 49 PetscCall(VecDuplicate(x,&y1)); 50 PetscCall(PetscObjectSetName((PetscObject) y1, "Frequency space vector")); 51 PetscCall(VecDuplicate(x,&y2)); 52 PetscCall(PetscObjectSetName((PetscObject) y2, "Frequency space window vector")); 53 PetscCall(VecDuplicate(x,&z1)); 54 PetscCall(PetscObjectSetName((PetscObject) z1, "Reconstructed convolution")); 55 PetscCall(VecDuplicate(x,&z2)); 56 PetscCall(PetscObjectSetName((PetscObject) z2, "Real space convolution")); 57 58 if (function == RANDOM) { 59 PetscCall(VecSetRandom(x, rdm)); 60 } else if (function == CONSTANT) { 61 PetscCall(VecSet(x, 1.0)); 62 } else if (function == TANH) { 63 PetscCall(VecGetArray(x, &a)); 64 for (i = 0; i < N; ++i) { 65 a[i] = tanh((i - N/2.0)*(10.0/N)); 66 } 67 PetscCall(VecRestoreArray(x, &a)); 68 } 69 if (view) PetscCall(VecView(x, PETSC_VIEWER_DRAW_WORLD)); 70 71 /* Create window function */ 72 PetscCall(VecGetArray(w, &a)); 73 for (i = 0; i < N; ++i) { 74 /* Step Function */ 75 a[i] = (i > N/4 && i < 3*N/4) ? 1.0 : 0.0; 76 /* Delta Function */ 77 /*a[i] = (i == N/2)? 1.0: 0.0; */ 78 } 79 PetscCall(VecRestoreArray(w, &a)); 80 if (view) PetscCall(VecView(w, PETSC_VIEWER_DRAW_WORLD)); 81 82 /* create FFTW object */ 83 PetscCall(MatCreateFFT(PETSC_COMM_SELF,DIM,dim,MATFFTW,&A)); 84 85 /* Convolve x with w*/ 86 PetscCall(MatMult(A,x,y1)); 87 PetscCall(MatMult(A,w,y2)); 88 PetscCall(VecPointwiseMult(y1, y1, y2)); 89 if (view && i == 0) PetscCall(VecView(y1, PETSC_VIEWER_DRAW_WORLD)); 90 PetscCall(MatMultTranspose(A,y1,z1)); 91 92 /* Compute the real space convolution */ 93 PetscCall(VecGetArray(x, &a)); 94 PetscCall(VecGetArray(w, &a2)); 95 PetscCall(VecGetArray(z2, &a3)); 96 for (i = 0; i < N; ++i) { 97 /* PetscInt checkInd = (i > N/2-1)? i-N/2: i+N/2;*/ 98 99 a3[i] = 0.0; 100 for (j = -N/2+1; j < N/2; ++j) { 101 PetscInt xpInd = (j < 0) ? N+j : j; 102 PetscInt diffInd = (i-j < 0) ? N-(j-i) : (i-j > N-1) ? i-j-N : i-j; 103 104 a3[i] += a[xpInd]*a2[diffInd]; 105 } 106 } 107 PetscCall(VecRestoreArray(x, &a)); 108 PetscCall(VecRestoreArray(w, &a2)); 109 PetscCall(VecRestoreArray(z2, &a3)); 110 111 /* compare z1 and z2. FFTW computes an unnormalized DFT, thus z1 = N*z2 */ 112 s = 1.0/(PetscReal)N; 113 PetscCall(VecScale(z1,s)); 114 if (view) PetscCall(VecView(z1, PETSC_VIEWER_DRAW_WORLD)); 115 if (view) PetscCall(VecView(z2, PETSC_VIEWER_DRAW_WORLD)); 116 PetscCall(VecAXPY(z1,-1.0,z2)); 117 PetscCall(VecNorm(z1,NORM_1,&enorm)); 118 if (enorm > 1.e-11) { 119 PetscCall(PetscPrintf(PETSC_COMM_SELF," Error norm of |z1 - z2| %g\n",(double)enorm)); 120 } 121 122 /* free spaces */ 123 PetscCall(VecDestroy(&x)); 124 PetscCall(VecDestroy(&y1)); 125 PetscCall(VecDestroy(&y2)); 126 PetscCall(VecDestroy(&z1)); 127 PetscCall(VecDestroy(&z2)); 128 PetscCall(VecDestroy(&w)); 129 PetscCall(MatDestroy(&A)); 130 } 131 PetscCall(PetscRandomDestroy(&rdm)); 132 PetscCall(PetscFinalize()); 133 return 0; 134 } 135 136 /*TEST 137 138 build: 139 requires: fftw complex 140 141 test: 142 output_file: output/ex121.out 143 TODO: Example or FFTW interface is broken 144 145 TEST*/ 146