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*d71ae5a4SJacob Faibussowitsch int main(int argc, char **args) 14*d71ae5a4SJacob Faibussowitsch { 15c4762a1bSJed Brown Mat A; /* FFT Matrix */ 16c4762a1bSJed Brown Vec x, y, z; /* Work vectors */ 17c4762a1bSJed Brown Vec x1, y1, z1; /* Duplicate vectors */ 18c4762a1bSJed Brown PetscInt i, k; /* for iterating over dimensions */ 19c4762a1bSJed Brown PetscRandom rdm; /* for creating random input */ 20c4762a1bSJed Brown PetscScalar a; /* used to scale output */ 21c4762a1bSJed Brown PetscReal enorm; /* norm for sanity check */ 22c4762a1bSJed Brown PetscInt n = 10, N = 1; /* FFT dimension params */ 23c4762a1bSJed Brown PetscInt DIM, dim[5]; /* FFT params */ 24c4762a1bSJed Brown 25327415f7SBarry Smith PetscFunctionBeginUser; 269566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 279566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 28c4762a1bSJed Brown 29c4762a1bSJed Brown /* To create random input vector */ 309566063dSJacob Faibussowitsch PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm)); 319566063dSJacob Faibussowitsch PetscCall(PetscRandomSetFromOptions(rdm)); 32c4762a1bSJed Brown 33c4762a1bSJed Brown /* Iterate over dimensions, use PETSc-FFTW interface */ 34c4762a1bSJed Brown for (i = 1; i < 5; i++) { 35c4762a1bSJed Brown DIM = i; 36c4762a1bSJed Brown N = 1; 379371c9d4SSatish Balay for (k = 0; k < i; k++) { 389371c9d4SSatish Balay dim[k] = n; 399371c9d4SSatish Balay N *= n; 409371c9d4SSatish Balay } 41c4762a1bSJed Brown 429566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n %" PetscInt_FMT " dimensions: FFTW on vector of size %" PetscInt_FMT " \n", DIM, N)); 43c4762a1bSJed Brown 44c4762a1bSJed Brown /* create FFTW object */ 459566063dSJacob Faibussowitsch PetscCall(MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A)); 46c4762a1bSJed Brown /* create vectors of length N */ 479566063dSJacob Faibussowitsch PetscCall(MatCreateVecsFFTW(A, &x, &y, &z)); 48c4762a1bSJed Brown 499566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector")); 509566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector")); 519566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector")); 52c4762a1bSJed Brown 53c4762a1bSJed Brown /* Test vector duplication*/ 549566063dSJacob Faibussowitsch PetscCall(VecDuplicate(x, &x1)); 559566063dSJacob Faibussowitsch PetscCall(VecDuplicate(y, &y1)); 569566063dSJacob Faibussowitsch PetscCall(VecDuplicate(z, &z1)); 57c4762a1bSJed Brown 58c4762a1bSJed Brown /* Set values of space vector x, copy to duplicate */ 599566063dSJacob Faibussowitsch PetscCall(VecSetRandom(x, rdm)); 609566063dSJacob Faibussowitsch PetscCall(VecCopy(x, x1)); 61c4762a1bSJed Brown 62c4762a1bSJed Brown /* Apply FFTW_FORWARD and FFTW_BACKWARD */ 639566063dSJacob Faibussowitsch PetscCall(MatMult(A, x, y)); 649566063dSJacob Faibussowitsch PetscCall(MatMultTranspose(A, y, z)); 65c4762a1bSJed Brown 66c4762a1bSJed Brown /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */ 679566063dSJacob Faibussowitsch PetscCall(MatMult(A, x1, y1)); 689566063dSJacob Faibussowitsch PetscCall(MatMultTranspose(A, y1, z1)); 69c4762a1bSJed Brown 70c4762a1bSJed Brown /* Compare x and z1. FFTW computes an unnormalized DFT, thus z1 = N*x */ 71c4762a1bSJed Brown a = 1.0 / (PetscReal)N; 729566063dSJacob Faibussowitsch PetscCall(VecScale(z1, a)); 739566063dSJacob Faibussowitsch PetscCall(VecAXPY(z1, -1.0, x)); 749566063dSJacob Faibussowitsch PetscCall(VecNorm(z1, NORM_1, &enorm)); 759566063dSJacob Faibussowitsch if (enorm > 1.e-9) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error norm of |x - z1| %g\n", enorm)); 76c4762a1bSJed Brown 77c4762a1bSJed Brown /* free spaces */ 789566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x1)); 799566063dSJacob Faibussowitsch PetscCall(VecDestroy(&y1)); 809566063dSJacob Faibussowitsch PetscCall(VecDestroy(&z1)); 819566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x)); 829566063dSJacob Faibussowitsch PetscCall(VecDestroy(&y)); 839566063dSJacob Faibussowitsch PetscCall(VecDestroy(&z)); 849566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 85c4762a1bSJed Brown } 86c4762a1bSJed Brown 879566063dSJacob Faibussowitsch PetscCall(PetscRandomDestroy(&rdm)); 889566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 89b122ec5aSJacob Faibussowitsch return 0; 90c4762a1bSJed Brown } 91c4762a1bSJed Brown 92c4762a1bSJed Brown /*TEST 93c4762a1bSJed Brown 94c4762a1bSJed Brown build: 95c4762a1bSJed Brown requires: fftw complex 96c4762a1bSJed Brown 97c4762a1bSJed Brown test: 98c4762a1bSJed Brown suffix: 2 99c4762a1bSJed Brown nsize : 4 100c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_ESTIMATE -n 16 101c4762a1bSJed Brown 102c4762a1bSJed Brown test: 103c4762a1bSJed Brown suffix: 3 104c4762a1bSJed Brown nsize : 2 105c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_MEASURE -n 12 106c4762a1bSJed Brown 107c4762a1bSJed Brown test: 108c4762a1bSJed Brown suffix: 4 109c4762a1bSJed Brown nsize : 2 110c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_PATIENT -n 10 111c4762a1bSJed Brown 112c4762a1bSJed Brown test: 113c4762a1bSJed Brown suffix: 5 114c4762a1bSJed Brown nsize : 1 115c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE -n 5 116c4762a1bSJed Brown 117c4762a1bSJed Brown TEST*/ 118