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