1 static char help[] = "Scalable algorithm for Connected Components problem.\n\ 2 Entails changing the MatMult() for this matrix.\n\n\n"; 3 4 #include <petscmat.h> 5 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> 9 10 /* 11 Paper with Ananth: Frbenius norm of band was good proxy, but really want to know the rank outside 12 13 LU for diagonal blocks must do shifting instead of pivoting, preferably shifting individual rows (like Pardiso) 14 15 Draw picture of flow of reordering 16 17 Measure Forbenius norm of the blocks being dropped by Truncated SPIKE (might be contaminated by pivoting in LU) 18 19 Report on using Florida matrices (Maxim, Murat) 20 */ 21 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 27 28 c = A * c 29 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 32 33 max_{components} depth of BFS tree for component 34 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 */ 40 41 /* Only isolated vertices get a 1 on the diagonal */ 42 PetscErrorCode CreateGraph(MPI_Comm comm, PetscInt testnum, Mat *A) 43 { 44 Mat G; 45 PetscErrorCode ierr; 46 47 PetscFunctionBegin; 48 ierr = MatCreate(comm, &G);CHKERRQ(ierr); 49 /* The identity matrix */ 50 switch (testnum) { 51 case 0: 52 { 53 Vec D; 54 55 ierr = MatSetSizes(G, PETSC_DETERMINE, PETSC_DETERMINE, 5, 5);CHKERRQ(ierr); 56 ierr = MatSetUp(G);CHKERRQ(ierr); 57 ierr = MatCreateVecs(G, &D, NULL);CHKERRQ(ierr); 58 ierr = VecSet(D, 1.0);CHKERRQ(ierr); 59 ierr = MatDiagonalSet(G, D, INSERT_VALUES);CHKERRQ(ierr); 60 ierr = VecDestroy(&D);CHKERRQ(ierr); 61 } 62 break; 63 case 1: 64 { 65 PetscScalar vals[3] = {1.0, 1.0, 1.0}; 66 PetscInt cols[3]; 67 PetscInt rStart, rEnd, row; 68 69 ierr = MatSetSizes(G, PETSC_DETERMINE, PETSC_DETERMINE, 5, 5);CHKERRQ(ierr); 70 ierr = MatSetFromOptions(G);CHKERRQ(ierr); 71 ierr = MatSeqAIJSetPreallocation(G, 2, NULL);CHKERRQ(ierr); 72 ierr = MatSetUp(G);CHKERRQ(ierr); 73 ierr = MatGetOwnershipRange(G, &rStart, &rEnd);CHKERRQ(ierr); 74 row = 0; 75 cols[0] = 0; cols[1] = 1; 76 if ((row >= rStart) && (row < rEnd)) {ierr = MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES);CHKERRQ(ierr);} 77 row = 1; 78 cols[0] = 0; cols[1] = 1; 79 if ((row >= rStart) && (row < rEnd)) {ierr = MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES);CHKERRQ(ierr);} 80 row = 2; 81 cols[0] = 2; cols[1] = 3; 82 if ((row >= rStart) && (row < rEnd)) {ierr = MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES);CHKERRQ(ierr);} 83 row = 3; 84 cols[0] = 3; cols[1] = 4; 85 if ((row >= rStart) && (row < rEnd)) {ierr = MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES);CHKERRQ(ierr);} 86 row = 4; 87 cols[0] = 4; cols[1] = 2; 88 if ((row >= rStart) && (row < rEnd)) {ierr = MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES);CHKERRQ(ierr);} 89 ierr = MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 90 ierr = MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 91 } 92 break; 93 case 2: 94 { 95 PetscScalar vals[3] = {1.0, 1.0, 1.0}; 96 PetscInt cols[3]; 97 PetscInt rStart, rEnd, row; 98 99 ierr = MatSetSizes(G, PETSC_DETERMINE, PETSC_DETERMINE, 5, 5);CHKERRQ(ierr); 100 ierr = MatSetFromOptions(G);CHKERRQ(ierr); 101 ierr = MatSeqAIJSetPreallocation(G, 2, NULL);CHKERRQ(ierr); 102 ierr = MatSetUp(G);CHKERRQ(ierr); 103 ierr = MatGetOwnershipRange(G, &rStart, &rEnd);CHKERRQ(ierr); 104 row = 0; 105 cols[0] = 0; cols[1] = 4; 106 if ((row >= rStart) && (row < rEnd)) {ierr = MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES);CHKERRQ(ierr);} 107 row = 1; 108 cols[0] = 1; cols[1] = 2; 109 if ((row >= rStart) && (row < rEnd)) {ierr = MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES);CHKERRQ(ierr);} 110 row = 2; 111 cols[0] = 2; cols[1] = 3; 112 if ((row >= rStart) && (row < rEnd)) {ierr = MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES);CHKERRQ(ierr);} 113 row = 3; 114 cols[0] = 3; cols[1] = 1; 115 if ((row >= rStart) && (row < rEnd)) {ierr = MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES);CHKERRQ(ierr);} 116 row = 4; 117 cols[0] = 0; cols[1] = 4; 118 if ((row >= rStart) && (row < rEnd)) {ierr = MatSetValues(G, 1, &row, 2, cols, vals, INSERT_VALUES);CHKERRQ(ierr);} 119 ierr = MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 120 ierr = MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 121 } 122 break; 123 default: 124 SETERRQ1(comm, PETSC_ERR_PLIB, "Unknown test %d", testnum); 125 } 126 *A = G; 127 PetscFunctionReturn(0); 128 } 129 130 int main(int argc, char **argv) 131 { 132 MPI_Comm comm; 133 Mat A; /* A graph */ 134 Vec c; /* A vector giving the component of each vertex */ 135 Vec cold; /* The vector c from the last iteration */ 136 PetscScalar *carray; 137 PetscInt testnum = 0; 138 PetscInt V, vStart, vEnd, v, n; 139 PetscMPIInt size; 140 PetscErrorCode ierr; 141 142 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 143 comm = PETSC_COMM_WORLD; 144 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 145 /* Use matrix to encode a graph */ 146 ierr = PetscOptionsGetInt(NULL,NULL, "-testnum", &testnum, NULL);CHKERRQ(ierr); 147 ierr = CreateGraph(comm, testnum, &A);CHKERRQ(ierr); 148 ierr = MatGetSize(A, &V, NULL);CHKERRQ(ierr); 149 /* Replace matrix-vector multiplication with one that calculates the minimum rather than the sum */ 150 if (size == 1) { 151 ierr = MatShellSetOperation(A, MATOP_MULT, (void (*)) MatMultMax_SeqAIJ);CHKERRQ(ierr); 152 } else { 153 Mat_MPIAIJ *a = (Mat_MPIAIJ *) A->data; 154 155 ierr = MatShellSetOperation(a->A, MATOP_MULT, (void (*)) MatMultMax_SeqAIJ);CHKERRQ(ierr); 156 ierr = MatShellSetOperation(a->B, MATOP_MULT, (void (*)) MatMultMax_SeqAIJ);CHKERRQ(ierr); 157 ierr = MatShellSetOperation(a->B, MATOP_MULT_ADD, (void (*)) MatMultAddMax_SeqAIJ);CHKERRQ(ierr); 158 } 159 /* Initialize each vertex as a separate component */ 160 ierr = MatCreateVecs(A, &c, NULL);CHKERRQ(ierr); 161 ierr = MatGetOwnershipRange(A, &vStart, &vEnd);CHKERRQ(ierr); 162 ierr = VecGetArray(c, &carray);CHKERRQ(ierr); 163 for (v = vStart; v < vEnd; ++v) { 164 carray[v-vStart] = v; 165 } 166 ierr = VecRestoreArray(c, &carray);CHKERRQ(ierr); 167 /* Preprocess in parallel to find local components */ 168 /* Multiply until c does not change */ 169 ierr = VecDuplicate(c, &cold);CHKERRQ(ierr); 170 for (v = 0; v < V; ++v) { 171 Vec cnew = cold; 172 PetscBool stop; 173 174 ierr = MatMult(A, c, cnew);CHKERRQ(ierr); 175 ierr = VecEqual(c, cnew, &stop);CHKERRQ(ierr); 176 if (stop) break; 177 cold = c; 178 c = cnew; 179 } 180 /* Report */ 181 ierr = VecUniqueEntries(c, &n, NULL);CHKERRQ(ierr); 182 ierr = PetscPrintf(comm, "Components: %d Iterations: %d\n", n, v);CHKERRQ(ierr); 183 ierr = VecView(c, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 184 /* Cleanup */ 185 ierr = VecDestroy(&c);CHKERRQ(ierr); 186 ierr = VecDestroy(&cold);CHKERRQ(ierr); 187 ierr = PetscFinalize(); 188 return ierr; 189 } 190