1: #include <petscdmda.h>
  2: #include <adolc/adalloc.h>

  4: /*
  5:    REQUIRES configuration of PETSc with option --download-adolc.

  7:    For documentation on ADOL-C, see
  8:      $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf
  9: */

 11: /*
 12:   Wrapper function for allocating contiguous memory in a 2d array

 14:   Input parameters:
 15:   m,n - number of rows and columns of array, respectively

 17:   Output parameter:
 18:   A   - pointer to array for which memory is allocated

 20:   Note: Only arrays of doubles are currently accounted for in ADOL-C's myalloc2 function.
 21: */
 22: template <class T>
 23: PetscErrorCode AdolcMalloc2(PetscInt m, PetscInt n, T **A[])
 24: {
 25:   PetscFunctionBegin;
 26:   *A = myalloc2(m, n);
 27:   PetscFunctionReturn(PETSC_SUCCESS);
 28: }

 30: /*
 31:   Wrapper function for freeing memory allocated with AdolcMalloc2

 33:   Input parameter:
 34:   A - array to free memory of

 36:   Note: Only arrays of doubles are currently accounted for in ADOL-C's myfree2 function.
 37: */
 38: template <class T>
 39: PetscErrorCode AdolcFree2(T **A)
 40: {
 41:   PetscFunctionBegin;
 42:   myfree2(A);
 43:   PetscFunctionReturn(PETSC_SUCCESS);
 44: }

 46: /*
 47:   Shift indices in an array of type T to endow it with ghost points.
 48:   (e.g. This works for arrays of adoubles or arrays (of structs) thereof.)

 50:   Input parameters:
 51:   da   - distributed array upon which variables are defined
 52:   cgs  - contiguously allocated 1-array with as many entries as there are
 53:          interior and ghost points, in total

 55:   Output parameter:
 56:   array - contiguously allocated array of the appropriate dimension with
 57:           ghost points, pointing to the 1-array
 58: */
 59: template <class T>
 60: PetscErrorCode GiveGhostPoints(DM da, T *cgs, void *array)
 61: {
 62:   PetscInt dim;

 64:   PetscFunctionBegin;
 65:   PetscCall(DMDAGetInfo(da, &dim, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
 66:   if (dim == 1) {
 67:     PetscCall(GiveGhostPoints1d(da, (T **)array));
 68:   } else if (dim == 2) {
 69:     PetscCall(GiveGhostPoints2d(da, cgs, (T ***)array));
 70:   } else PetscCheck(dim != 3, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "GiveGhostPoints3d not yet implemented"); // TODO
 71:   PetscFunctionReturn(PETSC_SUCCESS);
 72: }

 74: /*
 75:   Shift indices in a 1-array of type T to endow it with ghost points.
 76:   (e.g. This works for arrays of adoubles or arrays (of structs) thereof.)

 78:   Input parameters:
 79:   da  - distributed array upon which variables are defined

 81:   Output parameter:
 82:   a1d - contiguously allocated 1-array
 83: */
 84: template <class T>
 85: PetscErrorCode GiveGhostPoints1d(DM da, T *a1d[])
 86: {
 87:   PetscInt gxs;

 89:   PetscFunctionBegin;
 90:   PetscCall(DMDAGetGhostCorners(da, &gxs, NULL, NULL, NULL, NULL, NULL));
 91:   *a1d -= gxs;
 92:   PetscFunctionReturn(PETSC_SUCCESS);
 93: }

 95: /*
 96:   Shift indices in a 2-array of type T to endow it with ghost points.
 97:   (e.g. This works for arrays of adoubles or arrays (of structs) thereof.)

 99:   Input parameters:
100:   da  - distributed array upon which variables are defined
101:   cgs - contiguously allocated 1-array with as many entries as there are
102:         interior and ghost points, in total

104:   Output parameter:
105:   a2d - contiguously allocated 2-array with ghost points, pointing to the
106:         1-array
107: */
108: template <class T>
109: PetscErrorCode GiveGhostPoints2d(DM da, T *cgs, T **a2d[])
110: {
111:   PetscInt gxs, gys, gxm, gym;

113:   PetscFunctionBegin;
114:   PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, NULL, &gxm, &gym, NULL));
115:   for (PetscInt j = 0; j < gym; j++) (*a2d)[j] = cgs + j * gxm - gxs;
116:   *a2d -= gys;
117:   PetscFunctionReturn(PETSC_SUCCESS);
118: }

120: /*
121:   Create a rectangular sub-identity of the m x m identity matrix, as an array.

123:   Input parameters:
124:   n - number of (adjacent) rows to take in slice
125:   s - starting row index

127:   Output parameter:
128:   S - resulting n x m submatrix
129: */
130: template <class T>
131: PetscErrorCode Subidentity(PetscInt n, PetscInt s, T **S)
132: {
133:   PetscInt i;

135:   PetscFunctionBegin;
136:   for (i = 0; i < n; i++) S[i][i + s] = 1.;
137:   PetscFunctionReturn(PETSC_SUCCESS);
138: }

140: /*
141:   Create an identity matrix, as an array.

143:   Input parameter:
144:   n - number of rows/columns
145:   I - n x n array with memory pre-allocated
146: */
147: template <class T>
148: PetscErrorCode Identity(PetscInt n, T **I)
149: {
150:   PetscFunctionBegin;
151:   PetscCall(Subidentity(n, 0, I));
152:   PetscFunctionReturn(PETSC_SUCCESS);
153: }