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