Actual source code: bandwidth.c

      1: #include <petsc/private/matimpl.h>

      3: /*@
      4:   MatComputeBandwidth - Calculate the full bandwidth of the matrix, meaning the width 2k+1 where k diagonals on either side are sufficient to contain all the matrix nonzeros.

      6:   Collective

      8:   Input Parameters:
      9: + A        - The `Mat`
     10: - fraction - An optional percentage of the Frobenius norm of the matrix that the bandwidth should enclose

     12:   Output Parameter:
     13: . bw - The matrix bandwidth

     15:   Level: beginner

     17: .seealso: `DMPlexCreate()`, `DMPlexSetConeSize()`, `DMPlexSetChart()`
     18: @*/
     19: PetscErrorCode MatComputeBandwidth(Mat A, PetscReal fraction, PetscInt *bw)
     20: {
     21:   PetscInt lbw[2] = {0, 0}, gbw[2];
     22:   PetscInt rStart, rEnd, r;

     24:   PetscFunctionBegin;
     27:   PetscAssertPointer(bw, 3);
     28:   PetscCheck(!(fraction > 0.0) || !(fraction < 1.0), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "We do not yet support a fractional bandwidth");
     29:   PetscCall(MatGetOwnershipRange(A, &rStart, &rEnd));
     30:   for (r = rStart; r < rEnd; ++r) {
     31:     const PetscInt *cols;
     32:     PetscInt        ncols;

     34:     PetscCall(MatGetRow(A, r, &ncols, &cols, NULL));
     35:     if (ncols) {
     36:       lbw[0] = PetscMax(lbw[0], r - cols[0]);
     37:       lbw[1] = PetscMax(lbw[1], cols[ncols - 1] - r);
     38:     }
     39:     PetscCall(MatRestoreRow(A, r, &ncols, &cols, NULL));
     40:   }
     41:   PetscCallMPI(MPIU_Allreduce(lbw, gbw, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)A)));
     42:   *bw = 2 * PetscMax(gbw[0], gbw[1]) + 1;
     43:   PetscFunctionReturn(PETSC_SUCCESS);
     44: }