1: #include <petsc/private/matimpl.h>      /*I "petscmat.h"  I*/
  2: #include <petscsf.h>

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

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

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

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