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: }