1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
2
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.
5
6 Collective on Mat
7
8 Input Parameters:
9 + A - The Mat
10 - fraction - An optional percentage of the Frobenius norm of the matrix that the bandwidth should enclose
11
12 Output Parameter:
13 . bw - The matrix bandwidth
14
15 Level: beginner
16
17 .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
18 @*/
MatComputeBandwidth(Mat A,PetscReal fraction,PetscInt * bw)19 PetscErrorCode MatComputeBandwidth(Mat A, PetscReal fraction, PetscInt *bw)
20 {
21 PetscInt lbw[2] = {0, 0}, gbw[2];
22 PetscInt rStart, rEnd, r;
23 PetscErrorCode ierr;
24
25 PetscFunctionBegin;
26 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
27 PetscValidLogicalCollectiveReal(A,fraction,2);
28 PetscValidPointer(bw, 3);
29 if ((fraction > 0.0) && (fraction < 1.0)) SETERRQ(PetscObjectComm((PetscObject) A), PETSC_ERR_SUP, "We do not yet support a fractional bandwidth");
30 ierr = MatGetOwnershipRange(A, &rStart, &rEnd);CHKERRQ(ierr);
31 for (r = rStart; r < rEnd; ++r) {
32 const PetscInt *cols;
33 PetscInt ncols;
34
35 ierr = MatGetRow(A, r, &ncols, &cols, NULL);CHKERRQ(ierr);
36 if (ncols) {
37 lbw[0] = PetscMax(lbw[0], r - cols[0]);
38 lbw[1] = PetscMax(lbw[1], cols[ncols-1] - r);
39 }
40 ierr = MatRestoreRow(A, r, &ncols, &cols, NULL);CHKERRQ(ierr);
41 }
42 ierr = MPIU_Allreduce(lbw, gbw, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) A));CHKERRQ(ierr);
43 *bw = 2*PetscMax(gbw[0], gbw[1]) + 1;
44 PetscFunctionReturn(0);
45 }
46