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:   MatCreate(comm, &G);
 47:   /* The identity matrix */
 48:   switch (testnum) {
 49:   case 0: {
 50:     Vec D;

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

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

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

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

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

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

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