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

 97: PETSC_EXTERN PetscErrorCode MatISColoringTest(Mat A, ISColoring iscoloring)
 98: {
 99:   PetscInt        nn, c, i, j, M, N, nc, nnz, col, row;
100:   const PetscInt *cia, *cja, *cols;
101:   IS             *isis;
102:   MPI_Comm        comm;
103:   PetscMPIInt     size;
104:   PetscBool       done;
105:   PetscBT         table;

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

109:   PetscObjectGetComm((PetscObject)A, &comm);
110:   MPI_Comm_size(comm, &size);

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

116:   MatGetSize(A, &M, NULL);
117:   PetscBTCreate(M, &table);
118:   for (c = 0; c < nn; c++) { /* for each color */
119:     ISGetSize(isis[c], &nc);
120:     if (nc <= 1) continue;

122:     PetscBTMemzero(M, table);
123:     ISGetIndices(isis[c], &cols);
124:     for (j = 0; j < nc; j++) { /* for each column */
125:       col = cols[j];
126:       nnz = cia[col + 1] - cia[col];
127:       for (i = 0; i < nnz; i++) {
128:         row = cja[cia[col] + i];
130:       }
131:     }
132:     ISRestoreIndices(isis[c], &cols);
133:   }
134:   PetscBTDestroy(&table);

136:   MatRestoreColumnIJ(A, 1, PETSC_FALSE, PETSC_TRUE, NULL, &cia, &cja, &done);
137:   return 0;
138: }