Actual source code: sparse.cxx

      1: #include <petscmat.h>
      2: #include <adolc/adolc_sparse.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:   Basic printing for sparsity pattern

     14:   Input parameters:
     15:   comm     - MPI communicator
     16:   m        - number of rows

     18:   Output parameter:
     19:   sparsity - matrix sparsity pattern, typically computed using an ADOL-C function such as jac_pat
     20: */
     21: PetscErrorCode PrintSparsity(MPI_Comm comm, PetscInt m, unsigned int **sparsity)
     22: {
     23:   PetscFunctionBegin;
     24:   PetscCall(PetscPrintf(comm, "Sparsity pattern:\n"));
     25:   for (PetscInt i = 0; i < m; i++) {
     26:     PetscCall(PetscPrintf(comm, "\n %2d: ", i));
     27:     for (PetscInt j = 1; j <= (PetscInt)sparsity[i][0]; j++) PetscCall(PetscPrintf(comm, " %2d ", sparsity[i][j]));
     28:   }
     29:   PetscCall(PetscPrintf(comm, "\n\n"));
     30:   PetscFunctionReturn(PETSC_SUCCESS);
     31: }

     33: /*
     34:   Generate a seed matrix defining the partition of columns of a matrix by a particular coloring,
     35:   used for matrix compression

     37:   Input parameter:
     38:   iscoloring - the index set coloring to be used

     40:   Output parameter:
     41:   S          - the resulting seed matrix

     43:   Notes:
     44:   Before calling GenerateSeedMatrix, Seed should be allocated as a logically 2d array with number of
     45:   rows equal to the matrix to be compressed and number of columns equal to the number of colors used
     46:   in iscoloring.
     47: */
     48: PetscErrorCode GenerateSeedMatrix(ISColoring iscoloring, PetscScalar **S)
     49: {
     50:   IS             *is;
     51:   PetscInt        p, size;
     52:   const PetscInt *indices;

     54:   PetscFunctionBegin;
     55:   PetscCall(ISColoringGetIS(iscoloring, PETSC_USE_POINTER, &p, &is));
     56:   for (PetscInt colour = 0; colour < p; colour++) {
     57:     PetscCall(ISGetLocalSize(is[colour], &size));
     58:     PetscCall(ISGetIndices(is[colour], &indices));
     59:     for (PetscInt j = 0; j < size; j++) S[indices[j]][colour] = 1.;
     60:     PetscCall(ISRestoreIndices(is[colour], &indices));
     61:   }
     62:   PetscCall(ISColoringRestoreIS(iscoloring, PETSC_USE_POINTER, &is));
     63:   PetscFunctionReturn(PETSC_SUCCESS);
     64: }

     66: /*
     67:   Establish a look-up vector whose entries contain the colour used for that diagonal entry. Clearly
     68:   we require the number of dependent and independent variables to be equal in this case.
     69:   Input parameters:
     70:   S        - the seed matrix defining the coloring
     71:   sparsity - the sparsity pattern of the matrix to be recovered, typically computed using an ADOL-C
     72:              function, such as jac_pat or hess_pat
     73:   m        - the number of rows of Seed (and the matrix to be recovered)
     74:   p        - the number of colors used (also the number of columns in Seed)
     75:   Output parameter:
     76:   R        - the recovery vector to be used for de-compression
     77: */
     78: PetscErrorCode GenerateSeedMatrixPlusRecovery(ISColoring iscoloring, PetscScalar **S, PetscScalar *R)
     79: {
     80:   IS             *is;
     81:   PetscInt        p, size, colour, j;
     82:   const PetscInt *indices;

     84:   PetscFunctionBegin;
     85:   PetscCall(ISColoringGetIS(iscoloring, PETSC_USE_POINTER, &p, &is));
     86:   for (colour = 0; colour < p; colour++) {
     87:     PetscCall(ISGetLocalSize(is[colour], &size));
     88:     PetscCall(ISGetIndices(is[colour], &indices));
     89:     for (j = 0; j < size; j++) {
     90:       S[indices[j]][colour] = 1.;
     91:       R[indices[j]]         = colour;
     92:     }
     93:     PetscCall(ISRestoreIndices(is[colour], &indices));
     94:   }
     95:   PetscCall(ISColoringRestoreIS(iscoloring, PETSC_USE_POINTER, &is));
     96:   PetscFunctionReturn(PETSC_SUCCESS);
     97: }

     99: /*
    100:   Establish a look-up matrix whose entries contain the column coordinates of the corresponding entry
    101:   in a matrix which has been compressed using the coloring defined by some seed matrix

    103:   Input parameters:
    104:   S        - the seed matrix defining the coloring
    105:   sparsity - the sparsity pattern of the matrix to be recovered, typically computed using an ADOL-C
    106:              function, such as jac_pat or hess_pat
    107:   m        - the number of rows of Seed (and the matrix to be recovered)
    108:   p        - the number of colors used (also the number of columns in Seed)

    110:   Output parameter:
    111:   R        - the recovery matrix to be used for de-compression
    112: */
    113: PetscErrorCode GetRecoveryMatrix(PetscScalar **S, unsigned int **sparsity, PetscInt m, PetscInt p, PetscScalar **R)
    114: {
    115:   PetscInt i, j, k, colour;

    117:   PetscFunctionBegin;
    118:   for (i = 0; i < m; i++) {
    119:     for (colour = 0; colour < p; colour++) {
    120:       R[i][colour] = -1.;
    121:       for (k = 1; k <= (PetscInt)sparsity[i][0]; k++) {
    122:         j = (PetscInt)sparsity[i][k];
    123:         if (S[j][colour] == 1.) {
    124:           R[i][colour] = j;
    125:           break;
    126:         }
    127:       }
    128:     }
    129:   }
    130:   PetscFunctionReturn(PETSC_SUCCESS);
    131: }

    133: /*
    134:   Recover the values of a sparse matrix from a compressed format and insert these into a matrix

    136:   Input parameters:
    137:   mode - use INSERT_VALUES or ADD_VALUES, as required
    138:   m    - number of rows of matrix.
    139:   p    - number of colors used in the compression of J (also the number of columns of R)
    140:   R    - recovery matrix to use in the decompression procedure
    141:   C    - compressed matrix to recover values from
    142:   a    - shift value for implicit problems (select NULL or unity for explicit problems)

    144:   Output parameter:
    145:   A    - Mat to be populated with values from compressed matrix
    146: */
    147: PetscErrorCode RecoverJacobian(Mat A, InsertMode mode, PetscInt m, PetscInt p, PetscScalar **R, PetscScalar **C, PetscReal *a)
    148: {
    149:   PetscFunctionBegin;
    150:   for (PetscInt i = 0; i < m; i++) {
    151:     for (PetscInt colour = 0; colour < p; colour++) {
    152:       PetscInt j = (PetscInt)R[i][colour];
    153:       if (j != -1) {
    154:         if (a) C[i][colour] *= *a;
    155:         PetscCall(MatSetValues(A, 1, &i, 1, &j, &C[i][colour], mode));
    156:       }
    157:     }
    158:   }
    159:   PetscFunctionReturn(PETSC_SUCCESS);
    160: }

    162: /*
    163:   Recover the values of the local portion of a sparse matrix from a compressed format and insert
    164:   these into the local portion of a matrix

    166:   Input parameters:
    167:   mode - use INSERT_VALUES or ADD_VALUES, as required
    168:   m    - number of rows of matrix.
    169:   p    - number of colors used in the compression of J (also the number of columns of R)
    170:   R    - recovery matrix to use in the decompression procedure
    171:   C    - compressed matrix to recover values from
    172:   a    - shift value for implicit problems (select NULL or unity for explicit problems)

    174:   Output parameter:
    175:   A    - Mat to be populated with values from compressed matrix
    176: */
    177: PetscErrorCode RecoverJacobianLocal(Mat A, InsertMode mode, PetscInt m, PetscInt p, PetscScalar **R, PetscScalar **C, PetscReal *a)
    178: {
    179:   PetscFunctionBegin;
    180:   for (PetscInt i = 0; i < m; i++) {
    181:     for (PetscInt colour = 0; colour < p; colour++) {
    182:       PetscInt j = (PetscInt)R[i][colour];
    183:       if (j != -1) {
    184:         if (a) C[i][colour] *= *a;
    185:         PetscCall(MatSetValuesLocal(A, 1, &i, 1, &j, &C[i][colour], mode));
    186:       }
    187:     }
    188:   }
    189:   PetscFunctionReturn(PETSC_SUCCESS);
    190: }

    192: /*
    193:   Recover the diagonal of the Jacobian from its compressed matrix format

    195:   Input parameters:
    196:   mode - use INSERT_VALUES or ADD_VALUES, as required
    197:   m    - number of rows of matrix.
    198:   R    - recovery vector to use in the decompression procedure
    199:   C    - compressed matrix to recover values from
    200:   a    - shift value for implicit problems (select NULL or unity for explicit problems)

    202:   Output parameter:
    203:   diag - Vec to be populated with values from compressed matrix
    204: */
    205: PetscErrorCode RecoverDiagonal(Vec diag, InsertMode mode, PetscInt m, PetscScalar *R, PetscScalar **C, PetscReal *a)
    206: {
    207:   PetscFunctionBegin;
    208:   for (PetscInt i = 0; i < m; i++) {
    209:     PetscInt colour = (PetscInt)R[i];
    210:     if (a) C[i][colour] *= *a;
    211:     PetscCall(VecSetValues(diag, 1, &i, &C[i][colour], mode));
    212:   }
    213:   PetscFunctionReturn(PETSC_SUCCESS);
    214: }

    216: /*
    217:   Recover the local portion of the diagonal of the Jacobian from its compressed matrix format

    219:   Input parameters:
    220:   mode - use INSERT_VALUES or ADD_VALUES, as required
    221:   m    - number of rows of matrix.
    222:   R    - recovery vector to use in the decompression procedure
    223:   C    - compressed matrix to recover values from
    224:   a    - shift value for implicit problems (select NULL or unity for explicit problems)

    226:   Output parameter:
    227:   diag - Vec to be populated with values from compressed matrix
    228: */
    229: PetscErrorCode RecoverDiagonalLocal(Vec diag, InsertMode mode, PetscInt m, PetscScalar *R, PetscScalar **C, PetscReal *a)
    230: {
    231:   PetscFunctionBegin;
    232:   for (PetscInt i = 0; i < m; i++) {
    233:     PetscInt colour = (PetscInt)R[i];
    234:     if (a) C[i][colour] *= *a;
    235:     PetscCall(VecSetValuesLocal(diag, 1, &i, &C[i][colour], mode));
    236:   }
    237:   PetscFunctionReturn(PETSC_SUCCESS);
    238: }