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