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 PetscInt main(PetscInt 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 PetscErrorCode ierr; 26 27 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 28 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(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", "ex121", 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 /* create vectors of length N=n^DIM */ 44 for (i = 0, N = 1; i < DIM; i++) N *= dim[i]; 45 ierr = PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n",DIM,N);CHKERRQ(ierr); 46 ierr = VecCreateSeq(PETSC_COMM_SELF,N,&x);CHKERRQ(ierr); 47 ierr = PetscObjectSetName((PetscObject) x, "Real space vector");CHKERRQ(ierr); 48 ierr = VecDuplicate(x,&w);CHKERRQ(ierr); 49 ierr = PetscObjectSetName((PetscObject) w, "Window vector");CHKERRQ(ierr); 50 ierr = VecDuplicate(x,&y1);CHKERRQ(ierr); 51 ierr = PetscObjectSetName((PetscObject) y1, "Frequency space vector");CHKERRQ(ierr); 52 ierr = VecDuplicate(x,&y2);CHKERRQ(ierr); 53 ierr = PetscObjectSetName((PetscObject) y2, "Frequency space window vector");CHKERRQ(ierr); 54 ierr = VecDuplicate(x,&z1);CHKERRQ(ierr); 55 ierr = PetscObjectSetName((PetscObject) z1, "Reconstructed convolution");CHKERRQ(ierr); 56 ierr = VecDuplicate(x,&z2);CHKERRQ(ierr); 57 ierr = PetscObjectSetName((PetscObject) z2, "Real space convolution");CHKERRQ(ierr); 58 59 if (function == RANDOM) { 60 ierr = VecSetRandom(x, rdm);CHKERRQ(ierr); 61 } else if (function == CONSTANT) { 62 ierr = VecSet(x, 1.0);CHKERRQ(ierr); 63 } else if (function == TANH) { 64 ierr = VecGetArray(x, &a);CHKERRQ(ierr); 65 for (i = 0; i < N; ++i) { 66 a[i] = tanh((i - N/2.0)*(10.0/N)); 67 } 68 ierr = VecRestoreArray(x, &a);CHKERRQ(ierr); 69 } 70 if (view) {ierr = VecView(x, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);} 71 72 /* Create window function */ 73 ierr = VecGetArray(w, &a);CHKERRQ(ierr); 74 for (i = 0; i < N; ++i) { 75 /* Step Function */ 76 a[i] = (i > N/4 && i < 3*N/4) ? 1.0 : 0.0; 77 /* Delta Function */ 78 /*a[i] = (i == N/2)? 1.0: 0.0; */ 79 } 80 ierr = VecRestoreArray(w, &a);CHKERRQ(ierr); 81 if (view) {ierr = VecView(w, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);} 82 83 /* create FFTW object */ 84 ierr = MatCreateFFT(PETSC_COMM_SELF,DIM,dim,MATFFTW,&A);CHKERRQ(ierr); 85 86 /* Convolve x with w*/ 87 ierr = MatMult(A,x,y1);CHKERRQ(ierr); 88 ierr = MatMult(A,w,y2);CHKERRQ(ierr); 89 ierr = VecPointwiseMult(y1, y1, y2);CHKERRQ(ierr); 90 if (view && i == 0) {ierr = VecView(y1, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);} 91 ierr = MatMultTranspose(A,y1,z1);CHKERRQ(ierr); 92 93 /* Compute the real space convolution */ 94 ierr = VecGetArray(x, &a);CHKERRQ(ierr); 95 ierr = VecGetArray(w, &a2);CHKERRQ(ierr); 96 ierr = VecGetArray(z2, &a3);CHKERRQ(ierr); 97 for (i = 0; i < N; ++i) { 98 /* PetscInt checkInd = (i > N/2-1)? i-N/2: i+N/2;*/ 99 100 a3[i] = 0.0; 101 for (j = -N/2+1; j < N/2; ++j) { 102 PetscInt xpInd = (j < 0) ? N+j : j; 103 PetscInt diffInd = (i-j < 0) ? N-(j-i) : (i-j > N-1) ? i-j-N : i-j; 104 105 a3[i] += a[xpInd]*a2[diffInd]; 106 } 107 } 108 ierr = VecRestoreArray(x, &a);CHKERRQ(ierr); 109 ierr = VecRestoreArray(w, &a2);CHKERRQ(ierr); 110 ierr = VecRestoreArray(z2, &a3);CHKERRQ(ierr); 111 112 /* compare z1 and z2. FFTW computes an unnormalized DFT, thus z1 = N*z2 */ 113 s = 1.0/(PetscReal)N; 114 ierr = VecScale(z1,s);CHKERRQ(ierr); 115 if (view) {ierr = VecView(z1, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);} 116 if (view) {ierr = VecView(z2, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);} 117 ierr = VecAXPY(z1,-1.0,z2);CHKERRQ(ierr); 118 ierr = VecNorm(z1,NORM_1,&enorm);CHKERRQ(ierr); 119 if (enorm > 1.e-11) { 120 ierr = PetscPrintf(PETSC_COMM_SELF," Error norm of |z1 - z2| %g\n",(double)enorm);CHKERRQ(ierr); 121 } 122 123 /* free spaces */ 124 ierr = VecDestroy(&x);CHKERRQ(ierr); 125 ierr = VecDestroy(&y1);CHKERRQ(ierr); 126 ierr = VecDestroy(&y2);CHKERRQ(ierr); 127 ierr = VecDestroy(&z1);CHKERRQ(ierr); 128 ierr = VecDestroy(&z2);CHKERRQ(ierr); 129 ierr = VecDestroy(&w);CHKERRQ(ierr); 130 ierr = MatDestroy(&A);CHKERRQ(ierr); 131 } 132 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr); 133 ierr = PetscFinalize(); 134 return ierr; 135 } 136 137 138 /*TEST 139 140 build: 141 requires: fftw complex 142 143 test: 144 output_file: output/ex121.out 145 TODO: Example or FFTW interface is broken 146 147 TEST*/ 148