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