Actual source code: valid.c

  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: {
  8:   Mat            m=mc->mat;
  9:   PetscSF        etor,etoc;
 10:   PetscInt       s,e;
 11:   PetscInt       ncolors,nrows,ncols;
 12:   IS             *colors;
 13:   PetscInt       i,j,k,l;
 14:   PetscInt       *staterow,*statecol,*statespread;
 15:   PetscInt       nindices;
 16:   const PetscInt *indices;
 17:   PetscInt       dist=mc->dist;
 18:   const PetscInt *degrees;
 19:   PetscInt       *stateleafrow,*stateleafcol,nleafrows,nleafcols,idx,nentries,maxcolors;
 20:   MPI_Datatype   itype = MPIU_INT;

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

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

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

117:   ISColoringGetIS(iscoloring,PETSC_USE_POINTER,&nn,&isis);

119:   PetscObjectGetComm((PetscObject)A,&comm);
120:   MPI_Comm_size(comm,&size);

123:   MatGetColumnIJ(A,0,PETSC_FALSE,PETSC_FALSE,&N,&cia,&cja,&done);

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

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

148:   MatRestoreColumnIJ(A,1,PETSC_FALSE,PETSC_TRUE,NULL,&cia,&cja,&done);
149:   return 0;
150: }