Actual source code: ex170.c

  1: static char help[] = "Scalable algorithm for Connected Components problem.\n\
  2: Entails changing the MatMult() for this matrix.\n\n\n";

  4: #include <petscmat.h>

  6: PETSC_EXTERN PetscErrorCode MatMultMax_SeqAIJ(Mat, Vec, Vec);
  7: PETSC_EXTERN PetscErrorCode MatMultAddMax_SeqAIJ(Mat, Vec, Vec, Vec);
  8: #include <../src/mat/impls/aij/mpi/mpiaij.h>

 10: /*
 11:   Paper with Ananth: Frbenius norm of band was good proxy, but really want to know the rank outside

 13:   LU for diagonal blocks must do shifting instead of pivoting, preferably shifting individual rows (like Pardiso)

 15:   Draw picture of flow of reordering

 17:   Measure Forbenius norm of the blocks being dropped by Truncated SPIKE (might be contaminated by pivoting in LU)

 19:   Report on using Florida matrices (Maxim, Murat)
 20: */

 22: /*
 23: I have thought about how to do this. Here is a prototype algorithm. Let A be
 24: the adjacency matrix (0 or 1), and let each component be identified by the
 25: lowest numbered vertex in it. We initialize a vector c so that each vertex is
 26: a component, c_i = i. Now we act on c with A, using a special product

 28:   c = A * c

 30: where we replace addition with min. The fixed point of this operation is a vector
 31: c which is the component for each vertex. The number of iterates is

 33:   max_{components} depth of BFS tree for component

 35: We can accelerate this algorithm by preprocessing all locals domains using the
 36: same algorithm. Then the number of iterations is bounded the depth of the BFS
 37: tree for the graph on supervertices defined over local components, which is
 38: bounded by p. In practice, this should be very fast.
 39: */

 41: /* Only isolated vertices get a 1 on the diagonal */
 42: PetscErrorCode CreateGraph(MPI_Comm comm, PetscInt testnum, Mat *A)
 43: {
 44:   Mat G;

 46:   PetscFunctionBegin;
 47:   PetscCall(MatCreate(comm, &G));
 48:   /* The identity matrix */
 49:   switch (testnum) {
 50:   case 0: {
 51:     Vec D;

 53:     PetscCall(MatSetSizes(G, PETSC_DETERMINE, PETSC_DETERMINE, 5, 5));
 54:     PetscCall(MatSetUp(G));
 55:     PetscCall(MatCreateVecs(G, &D, NULL));
 56:     PetscCall(VecSet(D, 1.0));
 57:     PetscCall(MatDiagonalSet(G, D, INSERT_VALUES));
 58:     PetscCall(VecDestroy(&D));
 59:   } break;
 60:   case 1: {
 61:     PetscScalar vals[3] = {1.0, 1.0, 1.0};
 62:     PetscInt    cols[3];
 63:     PetscInt    rStart, rEnd, row;

 65:     PetscCall(MatSetSizes(G, PETSC_DETERMINE, PETSC_DETERMINE, 5, 5));
 66:     PetscCall(MatSetFromOptions(G));
 67:     PetscCall(MatSeqAIJSetPreallocation(G, 2, NULL));
 68:     PetscCall(MatSetUp(G));
 69:     PetscCall(MatGetOwnershipRange(G, &rStart, &rEnd));
 70:     row     = 0;
 71:     cols[0] = 0;
 72:     cols[1] = 1;
 73:     if ((row >= rStart) && (row < rEnd)) PetscCall(MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES));
 74:     row     = 1;
 75:     cols[0] = 0;
 76:     cols[1] = 1;
 77:     if ((row >= rStart) && (row < rEnd)) PetscCall(MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES));
 78:     row     = 2;
 79:     cols[0] = 2;
 80:     cols[1] = 3;
 81:     if ((row >= rStart) && (row < rEnd)) PetscCall(MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES));
 82:     row     = 3;
 83:     cols[0] = 3;
 84:     cols[1] = 4;
 85:     if ((row >= rStart) && (row < rEnd)) PetscCall(MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES));
 86:     row     = 4;
 87:     cols[0] = 4;
 88:     cols[1] = 2;
 89:     if ((row >= rStart) && (row < rEnd)) PetscCall(MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES));
 90:     PetscCall(MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY));
 91:     PetscCall(MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY));
 92:   } break;
 93:   case 2: {
 94:     PetscScalar vals[3] = {1.0, 1.0, 1.0};
 95:     PetscInt    cols[3];
 96:     PetscInt    rStart, rEnd, row;

 98:     PetscCall(MatSetSizes(G, PETSC_DETERMINE, PETSC_DETERMINE, 5, 5));
 99:     PetscCall(MatSetFromOptions(G));
100:     PetscCall(MatSeqAIJSetPreallocation(G, 2, NULL));
101:     PetscCall(MatSetUp(G));
102:     PetscCall(MatGetOwnershipRange(G, &rStart, &rEnd));
103:     row     = 0;
104:     cols[0] = 0;
105:     cols[1] = 4;
106:     if ((row >= rStart) && (row < rEnd)) PetscCall(MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES));
107:     row     = 1;
108:     cols[0] = 1;
109:     cols[1] = 2;
110:     if ((row >= rStart) && (row < rEnd)) PetscCall(MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES));
111:     row     = 2;
112:     cols[0] = 2;
113:     cols[1] = 3;
114:     if ((row >= rStart) && (row < rEnd)) PetscCall(MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES));
115:     row     = 3;
116:     cols[0] = 3;
117:     cols[1] = 1;
118:     if ((row >= rStart) && (row < rEnd)) PetscCall(MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES));
119:     row     = 4;
120:     cols[0] = 0;
121:     cols[1] = 4;
122:     if ((row >= rStart) && (row < rEnd)) PetscCall(MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES));
123:     PetscCall(MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY));
124:     PetscCall(MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY));
125:   } break;
126:   default:
127:     SETERRQ(comm, PETSC_ERR_PLIB, "Unknown test %d", testnum);
128:   }
129:   *A = G;
130:   PetscFunctionReturn(PETSC_SUCCESS);
131: }

133: int main(int argc, char **argv)
134: {
135:   MPI_Comm     comm;
136:   Mat          A;    /* A graph */
137:   Vec          c;    /* A vector giving the component of each vertex */
138:   Vec          cold; /* The vector c from the last iteration */
139:   PetscScalar *carray;
140:   PetscInt     testnum = 0;
141:   PetscInt     V, vStart, vEnd, v, n;
142:   PetscMPIInt  size;

144:   PetscFunctionBeginUser;
145:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
146:   comm = PETSC_COMM_WORLD;
147:   PetscCallMPI(MPI_Comm_size(comm, &size));
148:   /* Use matrix to encode a graph */
149:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-testnum", &testnum, NULL));
150:   PetscCall(CreateGraph(comm, testnum, &A));
151:   PetscCall(MatGetSize(A, &V, NULL));
152:   /* Replace matrix-vector multiplication with one that calculates the minimum rather than the sum */
153:   if (size == 1) {
154:     PetscCall(MatShellSetOperation(A, MATOP_MULT, (void(*))MatMultMax_SeqAIJ));
155:   } else {
156:     Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data;

158:     PetscCall(MatShellSetOperation(a->A, MATOP_MULT, (void(*))MatMultMax_SeqAIJ));
159:     PetscCall(MatShellSetOperation(a->B, MATOP_MULT, (void(*))MatMultMax_SeqAIJ));
160:     PetscCall(MatShellSetOperation(a->B, MATOP_MULT_ADD, (void(*))MatMultAddMax_SeqAIJ));
161:   }
162:   /* Initialize each vertex as a separate component */
163:   PetscCall(MatCreateVecs(A, &c, NULL));
164:   PetscCall(MatGetOwnershipRange(A, &vStart, &vEnd));
165:   PetscCall(VecGetArray(c, &carray));
166:   for (v = vStart; v < vEnd; ++v) carray[v - vStart] = v;
167:   PetscCall(VecRestoreArray(c, &carray));
168:   /* Preprocess in parallel to find local components */
169:   /* Multiply until c does not change */
170:   PetscCall(VecDuplicate(c, &cold));
171:   for (v = 0; v < V; ++v) {
172:     Vec       cnew = cold;
173:     PetscBool stop;

175:     PetscCall(MatMult(A, c, cnew));
176:     PetscCall(VecEqual(c, cnew, &stop));
177:     if (stop) break;
178:     cold = c;
179:     c    = cnew;
180:   }
181:   /* Report */
182:   PetscCall(VecUniqueEntries(c, &n, NULL));
183:   PetscCall(PetscPrintf(comm, "Components: %d Iterations: %d\n", n, v));
184:   PetscCall(VecView(c, PETSC_VIEWER_STDOUT_WORLD));
185:   /* Cleanup */
186:   PetscCall(VecDestroy(&c));
187:   PetscCall(VecDestroy(&cold));
188:   PetscCall(PetscFinalize());
189:   return 0;
190: }