1c4762a1bSJed Brown static char help[] = "Test duplication/destruction of FFTW vecs \n\n"; 2c4762a1bSJed Brown 3c4762a1bSJed Brown /* 4c4762a1bSJed Brown Compiling the code: 5c4762a1bSJed Brown This code uses the FFTW interface. 6c4762a1bSJed Brown Use one of the options below to configure: 7c4762a1bSJed Brown --with-fftw-dir=/.... or --download-fftw 8c4762a1bSJed Brown Usage: 9c4762a1bSJed Brown mpiexec -np <np> ./ex228 10c4762a1bSJed Brown */ 11c4762a1bSJed Brown 12c4762a1bSJed Brown #include <petscmat.h> 13*9371c9d4SSatish Balay int main(int argc, char **args) { 14c4762a1bSJed Brown Mat A; /* FFT Matrix */ 15c4762a1bSJed Brown Vec x, y, z; /* Work vectors */ 16c4762a1bSJed Brown Vec x1, y1, z1; /* Duplicate vectors */ 17c4762a1bSJed Brown PetscInt i, k; /* for iterating over dimensions */ 18c4762a1bSJed Brown PetscRandom rdm; /* for creating random input */ 19c4762a1bSJed Brown PetscScalar a; /* used to scale output */ 20c4762a1bSJed Brown PetscReal enorm; /* norm for sanity check */ 21c4762a1bSJed Brown PetscInt n = 10, N = 1; /* FFT dimension params */ 22c4762a1bSJed Brown PetscInt DIM, dim[5]; /* FFT params */ 23c4762a1bSJed Brown 24327415f7SBarry Smith PetscFunctionBeginUser; 259566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 269566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 27c4762a1bSJed Brown 28c4762a1bSJed Brown /* To create random input vector */ 299566063dSJacob Faibussowitsch PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm)); 309566063dSJacob Faibussowitsch PetscCall(PetscRandomSetFromOptions(rdm)); 31c4762a1bSJed Brown 32c4762a1bSJed Brown /* Iterate over dimensions, use PETSc-FFTW interface */ 33c4762a1bSJed Brown for (i = 1; i < 5; i++) { 34c4762a1bSJed Brown DIM = i; 35c4762a1bSJed Brown N = 1; 36*9371c9d4SSatish Balay for (k = 0; k < i; k++) { 37*9371c9d4SSatish Balay dim[k] = n; 38*9371c9d4SSatish Balay N *= n; 39*9371c9d4SSatish Balay } 40c4762a1bSJed Brown 419566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n %" PetscInt_FMT " dimensions: FFTW on vector of size %" PetscInt_FMT " \n", DIM, N)); 42c4762a1bSJed Brown 43c4762a1bSJed Brown /* create FFTW object */ 449566063dSJacob Faibussowitsch PetscCall(MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A)); 45c4762a1bSJed Brown /* create vectors of length N */ 469566063dSJacob Faibussowitsch PetscCall(MatCreateVecsFFTW(A, &x, &y, &z)); 47c4762a1bSJed Brown 489566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector")); 499566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector")); 509566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector")); 51c4762a1bSJed Brown 52c4762a1bSJed Brown /* Test vector duplication*/ 539566063dSJacob Faibussowitsch PetscCall(VecDuplicate(x, &x1)); 549566063dSJacob Faibussowitsch PetscCall(VecDuplicate(y, &y1)); 559566063dSJacob Faibussowitsch PetscCall(VecDuplicate(z, &z1)); 56c4762a1bSJed Brown 57c4762a1bSJed Brown /* Set values of space vector x, copy to duplicate */ 589566063dSJacob Faibussowitsch PetscCall(VecSetRandom(x, rdm)); 599566063dSJacob Faibussowitsch PetscCall(VecCopy(x, x1)); 60c4762a1bSJed Brown 61c4762a1bSJed Brown /* Apply FFTW_FORWARD and FFTW_BACKWARD */ 629566063dSJacob Faibussowitsch PetscCall(MatMult(A, x, y)); 639566063dSJacob Faibussowitsch PetscCall(MatMultTranspose(A, y, z)); 64c4762a1bSJed Brown 65c4762a1bSJed Brown /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */ 669566063dSJacob Faibussowitsch PetscCall(MatMult(A, x1, y1)); 679566063dSJacob Faibussowitsch PetscCall(MatMultTranspose(A, y1, z1)); 68c4762a1bSJed Brown 69c4762a1bSJed Brown /* Compare x and z1. FFTW computes an unnormalized DFT, thus z1 = N*x */ 70c4762a1bSJed Brown a = 1.0 / (PetscReal)N; 719566063dSJacob Faibussowitsch PetscCall(VecScale(z1, a)); 729566063dSJacob Faibussowitsch PetscCall(VecAXPY(z1, -1.0, x)); 739566063dSJacob Faibussowitsch PetscCall(VecNorm(z1, NORM_1, &enorm)); 749566063dSJacob Faibussowitsch if (enorm > 1.e-9) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error norm of |x - z1| %g\n", enorm)); 75c4762a1bSJed Brown 76c4762a1bSJed Brown /* free spaces */ 779566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x1)); 789566063dSJacob Faibussowitsch PetscCall(VecDestroy(&y1)); 799566063dSJacob Faibussowitsch PetscCall(VecDestroy(&z1)); 809566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x)); 819566063dSJacob Faibussowitsch PetscCall(VecDestroy(&y)); 829566063dSJacob Faibussowitsch PetscCall(VecDestroy(&z)); 839566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 84c4762a1bSJed Brown } 85c4762a1bSJed Brown 869566063dSJacob Faibussowitsch PetscCall(PetscRandomDestroy(&rdm)); 879566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 88b122ec5aSJacob Faibussowitsch return 0; 89c4762a1bSJed Brown } 90c4762a1bSJed Brown 91c4762a1bSJed Brown /*TEST 92c4762a1bSJed Brown 93c4762a1bSJed Brown build: 94c4762a1bSJed Brown requires: fftw complex 95c4762a1bSJed Brown 96c4762a1bSJed Brown test: 97c4762a1bSJed Brown suffix: 2 98c4762a1bSJed Brown nsize : 4 99c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_ESTIMATE -n 16 100c4762a1bSJed Brown 101c4762a1bSJed Brown test: 102c4762a1bSJed Brown suffix: 3 103c4762a1bSJed Brown nsize : 2 104c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_MEASURE -n 12 105c4762a1bSJed Brown 106c4762a1bSJed Brown test: 107c4762a1bSJed Brown suffix: 4 108c4762a1bSJed Brown nsize : 2 109c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_PATIENT -n 10 110c4762a1bSJed Brown 111c4762a1bSJed Brown test: 112c4762a1bSJed Brown suffix: 5 113c4762a1bSJed Brown nsize : 1 114c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE -n 5 115c4762a1bSJed Brown 116c4762a1bSJed Brown TEST*/ 117