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