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,°rees);
49: PetscSFComputeDegreeEnd(etor,°rees);
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,°rees);
70: PetscSFComputeDegreeEnd(etoc,°rees);
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: }