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