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