1:  #include <petsc/private/matimpl.h>
  2:  #include <petscsf.h>

  4: PETSC_EXTERN PetscErrorCode MatColoringCreateBipartiteGraph(MatColoring,PetscSF *,PetscSF *);

  6: PETSC_EXTERN PetscErrorCode MatColoringTestValid(MatColoring mc,ISColoring coloring)
  7: {
  9:   Mat            m=mc->mat;
 10:   PetscSF        etor,etoc;
 11:   PetscInt       s,e;
 12:   PetscInt       ncolors,nrows,ncols;
 13:   IS             *colors;
 14:   PetscInt       i,j,k,l;
 15:   PetscInt       *staterow,*statecol,*statespread;
 16:   PetscInt       nindices;
 17:   const PetscInt *indices;
 18:   PetscInt       dist=mc->dist;
 19:   const PetscInt *degrees;
 20:   PetscInt       *stateleafrow,*stateleafcol,nleafrows,nleafcols,idx,nentries,maxcolors;
 21:   MPI_Datatype   itype;

 24:   MatColoringGetMaxColors(mc,&maxcolors);
 25:   PetscDataTypeToMPIDataType(PETSC_INT,&itype);
 26:   /* get the communication structures and the colors */
 27:   MatColoringCreateBipartiteGraph(mc,&etoc,&etor);
 28:   ISColoringGetIS(coloring,&ncolors,&colors);
 29:   PetscSFGetGraph(etor,&nrows,&nleafrows,NULL,NULL);
 30:   PetscSFGetGraph(etoc,&ncols,&nleafcols,NULL,NULL);
 31:   MatGetOwnershipRangeColumn(m,&s,&e);
 32:   PetscMalloc1(ncols,&statecol);
 33:   PetscMalloc1(nrows,&staterow);
 34:   PetscMalloc1(nleafcols,&stateleafcol);
 35:   PetscMalloc1(nleafrows,&stateleafrow);

 37:   for (l=0;l<ncolors;l++) {
 38:     if (l > maxcolors) break;
 39:     for (k=0;k<ncols;k++) {
 40:       statecol[k] = -1;
 41:     }
 42:     ISGetLocalSize(colors[l],&nindices);
 43:     ISGetIndices(colors[l],&indices);
 44:     for (k=0;k<nindices;k++) {
 45:       statecol[indices[k]-s] = indices[k];
 46:     }
 47:     ISRestoreIndices(colors[l],&indices);
 48:     statespread = statecol;
 49:     for (k=0;k<dist;k++) {
 50:       if (k%2 == 1) {
 51:         PetscSFComputeDegreeBegin(etor,&degrees);
 52:         PetscSFComputeDegreeEnd(etor,&degrees);
 53:         nentries=0;
 54:         for(i=0;i<nrows;i++) {
 55:           nentries += degrees[i];
 56:         }
 57:         idx=0;
 58:         for(i=0;i<nrows;i++) {
 59:           for (j=0;j<degrees[i];j++) {
 60:             stateleafrow[idx] = staterow[i];
 61:             idx++;
 62:           }
 63:           statecol[i]=0.;
 64:         }
 65:         if (idx != nentries) SETERRQ2(PetscObjectComm((PetscObject)mc),PETSC_ERR_NOT_CONVERGED,"Bad number of entries %d vs %d",idx,nentries);
 66:         PetscLogEventBegin(MATCOLORING_Comm,mc,0,0,0);
 67:         PetscSFReduceBegin(etoc,itype,stateleafrow,statecol,MPI_MAX);
 68:         PetscSFReduceEnd(etoc,itype,stateleafrow,statecol,MPI_MAX);
 69:         PetscLogEventEnd(MATCOLORING_Comm,mc,0,0,0);
 70:         statespread = statecol;
 71:       } else {
 72:         PetscSFComputeDegreeBegin(etoc,&degrees);
 73:         PetscSFComputeDegreeEnd(etoc,&degrees);
 74:         nentries=0;
 75:         for(i=0;i<ncols;i++) {
 76:           nentries += degrees[i];
 77:         }
 78:         idx=0;
 79:         for(i=0;i<ncols;i++) {
 80:           for (j=0;j<degrees[i];j++) {
 81:             stateleafcol[idx] = statecol[i];
 82:             idx++;
 83:           }
 84:           staterow[i]=0.;
 85:         }
 86:         if (idx != nentries) SETERRQ2(PetscObjectComm((PetscObject)mc),PETSC_ERR_NOT_CONVERGED,"Bad number of entries %d vs %d",idx,nentries);
 87:         PetscLogEventBegin(MATCOLORING_Comm,mc,0,0,0);
 88:         PetscSFReduceBegin(etor,itype,stateleafcol,staterow,MPI_MAX);
 89:         PetscSFReduceEnd(etor,itype,stateleafcol,staterow,MPI_MAX);
 90:         PetscLogEventEnd(MATCOLORING_Comm,mc,0,0,0);
 91:         statespread = staterow;
 92:       }
 93:     }
 94:     for (k=0;k<nindices;k++) {
 95:       if (statespread[indices[k]-s] != indices[k]) {
 96:         PetscPrintf(PetscObjectComm((PetscObject)mc),"%d of color %d conflicts with %d\n",indices[k],l,statespread[indices[k]-s]);
 97:       }
 98:     }
 99:     ISRestoreIndices(colors[l],&indices);
100:   }
101:   PetscFree(statecol);
102:   PetscFree(staterow);
103:   PetscFree(stateleafcol);
104:   PetscFree(stateleafrow);
105:   PetscSFDestroy(&etor);
106:   PetscSFDestroy(&etoc);
107:   return(0);
108: }