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