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

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

  6: PETSC_EXTERN PetscErrorCode MatColoringTest(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 = MPIU_INT;

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

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

109: PETSC_EXTERN PetscErrorCode MatISColoringTest(Mat A,ISColoring iscoloring)
110: {
112:   PetscInt       nn,c,i,j,M,N,nc,nnz,col,row;
113:   const PetscInt *cia,*cja,*cols;
114:   IS             *isis;
115:   MPI_Comm       comm;
116:   PetscMPIInt    size;
117:   PetscBool      done;
118:   PetscBT        table;

121:   ISColoringGetIS(iscoloring,&nn,&isis);

123:   PetscObjectGetComm((PetscObject)A,&comm);
124:   MPI_Comm_size(comm,&size);
125:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only support sequential matrix");

127:   MatGetColumnIJ(A,0,PETSC_FALSE,PETSC_FALSE,&N,&cia,&cja,&done);
128:   if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");

130:   MatGetSize(A,&M,NULL);
131:   PetscBTCreate(M,&table);
132:   for (c=0; c<nn; c++) { /* for each color */
133:     ISGetSize(isis[c],&nc);
134:     if (nc <= 1) continue;

136:     PetscBTMemzero(M,table);
137:     ISGetIndices(isis[c],&cols);
138:     for (j=0; j<nc; j++) { /* for each column */
139:       col = cols[j];
140:       nnz = cia[col+1] - cia[col];
141:       for (i=0; i<nnz; i++) {
142:         row = cja[cia[col]+i];
143:         if (PetscBTLookupSet(table,row)) {
144:           SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"color %D, col %D: row %D already in this color",c,col,row);
145:         }
146:       }
147:     }
148:     ISRestoreIndices(isis[c],&cols);
149:   }
150:   PetscBTDestroy(&table);

152:   MatRestoreColumnIJ(A,1,PETSC_FALSE,PETSC_TRUE,NULL,&cia,&cja,&done);
153:   return(0);
154: }