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 @*/ 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