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 PetscInt dim; 60c4762a1bSJed Brown 61c4762a1bSJed Brown PetscFunctionBegin; 62*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0)); 63c4762a1bSJed Brown if (dim == 1) { 64*5f80ce2aSJacob Faibussowitsch CHKERRQ(GiveGhostPoints1d(da,(T**)array)); 65c4762a1bSJed Brown } else if (dim == 2) { 66*5f80ce2aSJacob Faibussowitsch CHKERRQ(GiveGhostPoints2d(da,cgs,(T***)array)); 672c71b3e2SJacob Faibussowitsch } else PetscCheckFalse(dim == 3,PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"GiveGhostPoints3d not yet implemented"); // TODO 68c4762a1bSJed Brown PetscFunctionReturn(0); 69c4762a1bSJed Brown } 70c4762a1bSJed Brown 71c4762a1bSJed Brown /* 72c4762a1bSJed Brown Shift indices in a 1-array of type T to endow it with ghost points. 73c4762a1bSJed Brown (e.g. This works for arrays of adoubles or arrays (of structs) thereof.) 74c4762a1bSJed Brown 75c4762a1bSJed Brown Input parameters: 76c4762a1bSJed Brown da - distributed array upon which variables are defined 77c4762a1bSJed Brown 78c4762a1bSJed Brown Output parameter: 79c4762a1bSJed Brown a1d - contiguously allocated 1-array 80c4762a1bSJed Brown */ 81c4762a1bSJed Brown template <class T> PetscErrorCode GiveGhostPoints1d(DM da,T *a1d[]) 82c4762a1bSJed Brown { 83c4762a1bSJed Brown PetscInt gxs; 84c4762a1bSJed Brown 85c4762a1bSJed Brown PetscFunctionBegin; 86*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDAGetGhostCorners(da,&gxs,NULL,NULL,NULL,NULL,NULL)); 87c4762a1bSJed Brown *a1d -= gxs; 88c4762a1bSJed Brown PetscFunctionReturn(0); 89c4762a1bSJed Brown } 90c4762a1bSJed Brown 91c4762a1bSJed Brown /* 92c4762a1bSJed Brown Shift indices in a 2-array of type T to endow it with ghost points. 93c4762a1bSJed Brown (e.g. This works for arrays of adoubles or arrays (of structs) thereof.) 94c4762a1bSJed Brown 95c4762a1bSJed Brown Input parameters: 96c4762a1bSJed Brown da - distributed array upon which variables are defined 97c4762a1bSJed Brown cgs - contiguously allocated 1-array with as many entries as there are 98c4762a1bSJed Brown interior and ghost points, in total 99c4762a1bSJed Brown 100c4762a1bSJed Brown Output parameter: 101c4762a1bSJed Brown a2d - contiguously allocated 2-array with ghost points, pointing to the 102c4762a1bSJed Brown 1-array 103c4762a1bSJed Brown */ 104c4762a1bSJed Brown template <class T> PetscErrorCode GiveGhostPoints2d(DM da,T *cgs,T **a2d[]) 105c4762a1bSJed Brown { 106*5f80ce2aSJacob Faibussowitsch PetscInt gxs,gys,gxm,gym; 107c4762a1bSJed Brown 108c4762a1bSJed Brown PetscFunctionBegin; 109*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMDAGetGhostCorners(da,&gxs,&gys,NULL,&gxm,&gym,NULL)); 110*5f80ce2aSJacob Faibussowitsch for (PetscInt j=0; j<gym; j++) (*a2d)[j] = cgs + j*gxm - gxs; 111c4762a1bSJed Brown *a2d -= gys; 112c4762a1bSJed Brown PetscFunctionReturn(0); 113c4762a1bSJed Brown } 114c4762a1bSJed Brown 115c4762a1bSJed Brown /* 116c4762a1bSJed Brown Create a rectangular sub-identity of the m x m identity matrix, as an array. 117c4762a1bSJed Brown 118c4762a1bSJed Brown Input parameters: 119c4762a1bSJed Brown n - number of (adjacent) rows to take in slice 120c4762a1bSJed Brown s - starting row index 121c4762a1bSJed Brown 122c4762a1bSJed Brown Output parameter: 123c4762a1bSJed Brown S - resulting n x m submatrix 124c4762a1bSJed Brown */ 125c4762a1bSJed Brown template <class T> PetscErrorCode Subidentity(PetscInt n,PetscInt s,T **S) 126c4762a1bSJed Brown { 127c4762a1bSJed Brown PetscInt i; 128c4762a1bSJed Brown 129c4762a1bSJed Brown PetscFunctionBegin; 130c4762a1bSJed Brown for (i=0; i<n; i++) { 131c4762a1bSJed Brown S[i][i+s] = 1.; 132c4762a1bSJed Brown } 133c4762a1bSJed Brown PetscFunctionReturn(0); 134c4762a1bSJed Brown } 135c4762a1bSJed Brown 136c4762a1bSJed Brown /* 137c4762a1bSJed Brown Create an identity matrix, as an array. 138c4762a1bSJed Brown 139c4762a1bSJed Brown Input parameter: 140c4762a1bSJed Brown n - number of rows/columns 141c4762a1bSJed Brown I - n x n array with memory pre-allocated 142c4762a1bSJed Brown */ 143c4762a1bSJed Brown template <class T> PetscErrorCode Identity(PetscInt n,T **I) 144c4762a1bSJed Brown { 145c4762a1bSJed Brown PetscFunctionBegin; 146*5f80ce2aSJacob Faibussowitsch CHKERRQ(Subidentity(n,0,I)); 147c4762a1bSJed Brown PetscFunctionReturn(0); 148c4762a1bSJed Brown } 149