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 */ 22*9371c9d4SSatish Balay template <class T> 23*9371c9d4SSatish Balay PetscErrorCode AdolcMalloc2(PetscInt m, PetscInt n, T **A[]) { 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 */ 37*9371c9d4SSatish Balay template <class T> 38*9371c9d4SSatish Balay PetscErrorCode AdolcFree2(T **A) { 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 */ 57*9371c9d4SSatish Balay template <class T> 58*9371c9d4SSatish Balay PetscErrorCode GiveGhostPoints(DM da, T *cgs, void *array) { 59c4762a1bSJed Brown PetscInt dim; 60c4762a1bSJed Brown 61c4762a1bSJed Brown PetscFunctionBegin; 629566063dSJacob Faibussowitsch PetscCall(DMDAGetInfo(da, &dim, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 63c4762a1bSJed Brown if (dim == 1) { 649566063dSJacob Faibussowitsch PetscCall(GiveGhostPoints1d(da, (T **)array)); 65c4762a1bSJed Brown } else if (dim == 2) { 669566063dSJacob Faibussowitsch PetscCall(GiveGhostPoints2d(da, cgs, (T ***)array)); 6708401ef6SPierre Jolivet } else PetscCheck(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 */ 81*9371c9d4SSatish Balay template <class T> 82*9371c9d4SSatish Balay PetscErrorCode GiveGhostPoints1d(DM da, T *a1d[]) { 83c4762a1bSJed Brown PetscInt gxs; 84c4762a1bSJed Brown 85c4762a1bSJed Brown PetscFunctionBegin; 869566063dSJacob Faibussowitsch PetscCall(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 */ 104*9371c9d4SSatish Balay template <class T> 105*9371c9d4SSatish Balay PetscErrorCode GiveGhostPoints2d(DM da, T *cgs, T **a2d[]) { 1065f80ce2aSJacob Faibussowitsch PetscInt gxs, gys, gxm, gym; 107c4762a1bSJed Brown 108c4762a1bSJed Brown PetscFunctionBegin; 1099566063dSJacob Faibussowitsch PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, NULL, &gxm, &gym, NULL)); 1105f80ce2aSJacob 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 */ 125*9371c9d4SSatish Balay template <class T> 126*9371c9d4SSatish Balay PetscErrorCode Subidentity(PetscInt n, PetscInt s, T **S) { 127c4762a1bSJed Brown PetscInt i; 128c4762a1bSJed Brown 129c4762a1bSJed Brown PetscFunctionBegin; 130*9371c9d4SSatish Balay for (i = 0; i < n; i++) { S[i][i + s] = 1.; } 131c4762a1bSJed Brown PetscFunctionReturn(0); 132c4762a1bSJed Brown } 133c4762a1bSJed Brown 134c4762a1bSJed Brown /* 135c4762a1bSJed Brown Create an identity matrix, as an array. 136c4762a1bSJed Brown 137c4762a1bSJed Brown Input parameter: 138c4762a1bSJed Brown n - number of rows/columns 139c4762a1bSJed Brown I - n x n array with memory pre-allocated 140c4762a1bSJed Brown */ 141*9371c9d4SSatish Balay template <class T> 142*9371c9d4SSatish Balay PetscErrorCode Identity(PetscInt n, T **I) { 143c4762a1bSJed Brown PetscFunctionBegin; 1449566063dSJacob Faibussowitsch PetscCall(Subidentity(n, 0, I)); 145c4762a1bSJed Brown PetscFunctionReturn(0); 146c4762a1bSJed Brown } 147