[petsc-dev] [petsc-maint] Feature: allow unsorted row IS in MatGetSubMatrix_SeqAIJ()

[Barry Smith]

   Should actually send these things to pets-dev

   The docs for MatGetSubMatrix() have the line

"The rows in isrow will be sorted into the same order as the original matrix on each process."

    The docs for MatGetSubMatrices() say 

   "Currently both row and column indices must be sorted to guarantee
   correctness with all matrix types."

PetscErrorCode MatGetSubMatrices_MPIAIJ(Mat C,PetscInt ismax,const IS isrow[],const IS iscol[],MatReuse scall,Mat *submat[])
{
  PetscErrorCode ierr;
  PetscInt       nmax,nstages_local,nstages,i,pos,max_no,nrow,ncol;
  PetscBool      rowflag,colflag,wantallmatrix=PETSC_FALSE,twantallmatrix,*allcolumns;

  PetscFunctionBegin;
  /* Currently, unsorted column indices will result in inverted column indices in the resulting submatrices. */
  /* It would make sense to error out in MatGetSubMatrices_MPIAIJ_Local(), the most impl-specific level.
     However, there are more careful users of MatGetSubMatrices_MPIAIJ_Local() -- MatPermute_MPIAIJ() -- that
     take care to order the result correctly by assembling it with MatSetValues() (after preallocating).
   */
  for (i = 0; i < ismax; ++i) {
    PetscBool sorted;
    ierr = ISSorted(iscol[i], &sorted);CHKERRQ(ierr);
    if (!sorted) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Column index set %D not sorted", i);
  }

  In other words the sorted business seems a little wacky at the moment.  I'm not sure what to conclude.


   Barry

> On Jan 19, 2015, at 11:36 AM, Tobin Isaac <tisaac at ices.utexas.edu> wrote:
> 
> 
> There doesn't seem to be any reason for MatGetSubMatrix_SeqAIJ() to
> require the row IS to be sorted, and this requirement differs from
> other implementations, e.g. _MPIAIJ.  I have a feature branch that
> removes the requirement.  Can I merge it into next for testing?
> 
> Thanks,
>  Toby