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