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