1c4762a1bSJed Brown #include <petscdmda.h> 2c4762a1bSJed Brown #include <adolc/adalloc.h> 3c4762a1bSJed Brown 4c4762a1bSJed Brown /* 5c4762a1bSJed Brown REQUIRES configuration of PETSc with option --download-adolc. 6c4762a1bSJed Brown 7c4762a1bSJed Brown For documentation on ADOL-C, see 8c4762a1bSJed Brown $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf 9c4762a1bSJed Brown */ 10c4762a1bSJed Brown 11c4762a1bSJed Brown /* 12c4762a1bSJed Brown Wrapper function for allocating contiguous memory in a 2d array 13c4762a1bSJed Brown 14c4762a1bSJed Brown Input parameters: 15c4762a1bSJed Brown m,n - number of rows and columns of array, respectively 16c4762a1bSJed Brown 17c4762a1bSJed Brown Output parameter: 18c4762a1bSJed Brown A - pointer to array for which memory is allocated 19c4762a1bSJed Brown 20c4762a1bSJed Brown Note: Only arrays of doubles are currently accounted for in ADOL-C's myalloc2 function. 21c4762a1bSJed Brown */ 22c4762a1bSJed Brown template <class T> PetscErrorCode AdolcMalloc2(PetscInt m,PetscInt n,T **A[]) 23c4762a1bSJed Brown { 24c4762a1bSJed Brown PetscFunctionBegin; 25c4762a1bSJed Brown *A = myalloc2(m,n); 26c4762a1bSJed Brown PetscFunctionReturn(0); 27c4762a1bSJed Brown } 28c4762a1bSJed Brown 29c4762a1bSJed Brown /* 30c4762a1bSJed Brown Wrapper function for freeing memory allocated with AdolcMalloc2 31c4762a1bSJed Brown 32c4762a1bSJed Brown Input parameter: 33c4762a1bSJed Brown A - array to free memory of 34c4762a1bSJed Brown 35c4762a1bSJed Brown Note: Only arrays of doubles are currently accounted for in ADOL-C's myfree2 function. 36c4762a1bSJed Brown */ 37c4762a1bSJed Brown template <class T> PetscErrorCode AdolcFree2(T **A) 38c4762a1bSJed Brown { 39c4762a1bSJed Brown PetscFunctionBegin; 40c4762a1bSJed Brown myfree2(A); 41c4762a1bSJed Brown PetscFunctionReturn(0); 42c4762a1bSJed Brown } 43c4762a1bSJed Brown 44c4762a1bSJed Brown /* 45c4762a1bSJed Brown Shift indices in an array of type T to endow it with ghost points. 46c4762a1bSJed Brown (e.g. This works for arrays of adoubles or arrays (of structs) thereof.) 47c4762a1bSJed Brown 48c4762a1bSJed Brown Input parameters: 49c4762a1bSJed Brown da - distributed array upon which variables are defined 50c4762a1bSJed Brown cgs - contiguously allocated 1-array with as many entries as there are 51c4762a1bSJed Brown interior and ghost points, in total 52c4762a1bSJed Brown 53c4762a1bSJed Brown Output parameter: 54c4762a1bSJed Brown array - contiguously allocated array of the appropriate dimension with 55c4762a1bSJed Brown ghost points, pointing to the 1-array 56c4762a1bSJed Brown */ 57c4762a1bSJed Brown template <class T> PetscErrorCode GiveGhostPoints(DM da,T *cgs,void *array) 58c4762a1bSJed Brown { 59c4762a1bSJed Brown PetscErrorCode ierr; 60c4762a1bSJed Brown PetscInt dim; 61c4762a1bSJed Brown 62c4762a1bSJed Brown PetscFunctionBegin; 63c4762a1bSJed Brown ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 64c4762a1bSJed Brown if (dim == 1) { 65c4762a1bSJed Brown ierr = GiveGhostPoints1d(da,(T**)array);CHKERRQ(ierr); 66c4762a1bSJed Brown } else if (dim == 2) { 67c4762a1bSJed Brown ierr = GiveGhostPoints2d(da,cgs,(T***)array);CHKERRQ(ierr); 68*2c71b3e2SJacob Faibussowitsch } else PetscCheckFalse(dim == 3,PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"GiveGhostPoints3d not yet implemented"); // TODO 69c4762a1bSJed Brown PetscFunctionReturn(0); 70c4762a1bSJed Brown } 71c4762a1bSJed Brown 72c4762a1bSJed Brown /* 73c4762a1bSJed Brown Shift indices in a 1-array of type T to endow it with ghost points. 74c4762a1bSJed Brown (e.g. This works for arrays of adoubles or arrays (of structs) thereof.) 75c4762a1bSJed Brown 76c4762a1bSJed Brown Input parameters: 77c4762a1bSJed Brown da - distributed array upon which variables are defined 78c4762a1bSJed Brown 79c4762a1bSJed Brown Output parameter: 80c4762a1bSJed Brown a1d - contiguously allocated 1-array 81c4762a1bSJed Brown */ 82c4762a1bSJed Brown template <class T> PetscErrorCode GiveGhostPoints1d(DM da,T *a1d[]) 83c4762a1bSJed Brown { 84c4762a1bSJed Brown PetscErrorCode ierr; 85c4762a1bSJed Brown PetscInt gxs; 86c4762a1bSJed Brown 87c4762a1bSJed Brown PetscFunctionBegin; 88c4762a1bSJed Brown ierr = DMDAGetGhostCorners(da,&gxs,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 89c4762a1bSJed Brown *a1d -= gxs; 90c4762a1bSJed Brown PetscFunctionReturn(0); 91c4762a1bSJed Brown } 92c4762a1bSJed Brown 93c4762a1bSJed Brown /* 94c4762a1bSJed Brown Shift indices in a 2-array of type T to endow it with ghost points. 95c4762a1bSJed Brown (e.g. This works for arrays of adoubles or arrays (of structs) thereof.) 96c4762a1bSJed Brown 97c4762a1bSJed Brown Input parameters: 98c4762a1bSJed Brown da - distributed array upon which variables are defined 99c4762a1bSJed Brown cgs - contiguously allocated 1-array with as many entries as there are 100c4762a1bSJed Brown interior and ghost points, in total 101c4762a1bSJed Brown 102c4762a1bSJed Brown Output parameter: 103c4762a1bSJed Brown a2d - contiguously allocated 2-array with ghost points, pointing to the 104c4762a1bSJed Brown 1-array 105c4762a1bSJed Brown */ 106c4762a1bSJed Brown template <class T> PetscErrorCode GiveGhostPoints2d(DM da,T *cgs,T **a2d[]) 107c4762a1bSJed Brown { 108c4762a1bSJed Brown PetscErrorCode ierr; 109c4762a1bSJed Brown PetscInt gxs,gys,gxm,gym,j; 110c4762a1bSJed Brown 111c4762a1bSJed Brown PetscFunctionBegin; 112c4762a1bSJed Brown ierr = DMDAGetGhostCorners(da,&gxs,&gys,NULL,&gxm,&gym,NULL);CHKERRQ(ierr); 113c4762a1bSJed Brown for (j=0; j<gym; j++) 114c4762a1bSJed Brown (*a2d)[j] = cgs + j*gxm - gxs; 115c4762a1bSJed Brown *a2d -= gys; 116c4762a1bSJed Brown PetscFunctionReturn(0); 117c4762a1bSJed Brown } 118c4762a1bSJed Brown 119c4762a1bSJed Brown /* 120c4762a1bSJed Brown Create a rectangular sub-identity of the m x m identity matrix, as an array. 121c4762a1bSJed Brown 122c4762a1bSJed Brown Input parameters: 123c4762a1bSJed Brown n - number of (adjacent) rows to take in slice 124c4762a1bSJed Brown s - starting row index 125c4762a1bSJed Brown 126c4762a1bSJed Brown Output parameter: 127c4762a1bSJed Brown S - resulting n x m submatrix 128c4762a1bSJed Brown */ 129c4762a1bSJed Brown template <class T> PetscErrorCode Subidentity(PetscInt n,PetscInt s,T **S) 130c4762a1bSJed Brown { 131c4762a1bSJed Brown PetscInt i; 132c4762a1bSJed Brown 133c4762a1bSJed Brown PetscFunctionBegin; 134c4762a1bSJed Brown for (i=0; i<n; i++) { 135c4762a1bSJed Brown S[i][i+s] = 1.; 136c4762a1bSJed Brown } 137c4762a1bSJed Brown PetscFunctionReturn(0); 138c4762a1bSJed Brown } 139c4762a1bSJed Brown 140c4762a1bSJed Brown /* 141c4762a1bSJed Brown Create an identity matrix, as an array. 142c4762a1bSJed Brown 143c4762a1bSJed Brown Input parameter: 144c4762a1bSJed Brown n - number of rows/columns 145c4762a1bSJed Brown I - n x n array with memory pre-allocated 146c4762a1bSJed Brown */ 147c4762a1bSJed Brown template <class T> PetscErrorCode Identity(PetscInt n,T **I) 148c4762a1bSJed Brown { 149c4762a1bSJed Brown PetscErrorCode ierr; 150c4762a1bSJed Brown 151c4762a1bSJed Brown PetscFunctionBegin; 152c4762a1bSJed Brown ierr = Subidentity(n,0,I);CHKERRQ(ierr); 153c4762a1bSJed Brown PetscFunctionReturn(0); 154c4762a1bSJed Brown } 155