Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_Incidence_Matrix: Given the adjacency matrix of an undirected graph with no
3 : // self-loops, builds its corresponding incidence matrix
4 : //------------------------------------------------------------------------------
5 :
6 : // SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
7 : // SPDX-License-Identifier: Apache-2.0
8 : // See additional acknowledgments in the LICENSE file,
9 : // or contact permission@sei.cmu.edu for the full terms.
10 :
11 : // Contributed by Vidith Madhu, Texas A&M University
12 :
13 : // TODO: not ready for src; should handle all builtin types, with option to
14 : // typecast to INT64, UINT64, or FP64 as is currently done.
15 :
16 : // TODO: this method is required for MaximalMatching and CoarsenMatching
17 :
18 : // Given an undirected graph G, construct the incidence matrix E.
19 : /*
20 : The incidence matrix E has size n-by-e where the
21 : undirected graph G has n nodes and e edges. If the kth edge of G is the edge
22 : (i,j), then the column E(:,k) contains two entries: E(i,k) and E(j,k), which
23 : have the same value. If the graph G is weighted, then both E(i,k) and E(j,k)
24 : are equal to the weight of the (i,j) edge. If G is unweighted, then both are
25 : equal to 1 (and the matrix E is thus iso-valued).
26 :
27 : The type of the result is always compatible with the type of the input graph, but
28 : may be of a larger size; UINT8, UINT16, UINT32, INT8, INT16, INT32, INT64, and BOOL becomes
29 : a UINT64. UINT64 remains as a UINT64. FP32 and FP64 both become FP64.
30 :
31 : Note that complex types are NOT supported.
32 : */
33 :
34 : #include "LG_internal.h"
35 : #include "LAGraphX.h"
36 :
37 : #define LOADTRICKIM
38 : //#define dbg
39 :
40 : #undef LG_FREE_ALL
41 : #define LG_FREE_ALL \
42 : { \
43 : LAGraph_Free ((void**)(&row_indices), msg) ; \
44 : LAGraph_Free ((void**)(&col_indices), msg) ; \
45 : LAGraph_Free ((void**)(&values), msg) ; \
46 : LAGraph_Free ((void**)(&ramp), msg) ; \
47 : GrB_free (&con) ; \
48 : GrB_free (&E_half) ; \
49 : GrB_free (&A_tril) ; \
50 : GrB_free (&x); \
51 : GrB_free (&i); \
52 : GrB_free (&j); \
53 : GrB_free (&Ep); \
54 : GrB_free (&Ei); \
55 : GrB_free (&Ex); \
56 : GrB_free (&fullx); \
57 : GrB_free(&build_desc); \
58 : }
59 :
60 :
61 66 : int LAGraph_Incidence_Matrix
62 : (
63 : GrB_Matrix *result, // incidence
64 : LAGraph_Graph G, // must be undirected, no self-loops
65 : char *msg
66 : )
67 : {
68 : #if LAGRAPH_SUITESPARSE
69 :
70 66 : GrB_Matrix E = NULL ;
71 66 : GrB_Matrix E_half = NULL ;
72 66 : GrB_Matrix A_tril = NULL ;
73 66 : GrB_Vector i = NULL, j = NULL, x = NULL, fullx = NULL;
74 66 : GrB_Vector Ep = NULL, Ei = NULL, Ex = NULL;
75 66 : GrB_Descriptor build_desc = NULL;
76 66 : GrB_Index *row_indices = NULL ;
77 66 : GrB_Index *col_indices = NULL ;
78 66 : void *values = NULL ;
79 : #if GxB_IMPLEMENTATION < GxB_VERSION (10,0,0)
80 : GrB_Vector con = NULL; // unused, except by LG_FREE_ALL above
81 : #else
82 66 : GxB_Container con = NULL;
83 : #endif
84 :
85 66 : GrB_Index *ramp = NULL ;
86 :
87 66 : LG_ASSERT_MSG (
88 : G->kind == LAGraph_ADJACENCY_UNDIRECTED,
89 : LAGRAPH_INVALID_GRAPH,
90 : "G must be undirected"
91 : ) ;
92 :
93 66 : LG_ASSERT_MSG (G->nself_edges == 0, LAGRAPH_NO_SELF_EDGES_ALLOWED, "G->nself_edges must be zero") ;
94 :
95 66 : GrB_Matrix A = G->A ;
96 :
97 : char typename[LAGRAPH_MAX_NAME_LEN] ;
98 : GrB_Type type ;
99 : #if GxB_IMPLEMENTATION < GxB_VERSION (10,0,0)
100 : LG_TRY (LAGraph_Matrix_TypeName (typename, A, msg)) ;
101 : LG_TRY (LAGraph_TypeFromName (&type, typename, msg)) ;
102 : #else
103 : // No longer historical
104 66 : GRB_TRY (GxB_Matrix_type(&type, A));
105 : #endif
106 :
107 : GrB_Index nvals ;
108 : GrB_Index num_nodes, num_edges ;
109 :
110 66 : GRB_TRY (GrB_Matrix_nvals (&nvals, A)) ;
111 66 : num_edges = nvals / 2 ;
112 66 : GRB_TRY (GrB_Matrix_nrows (&num_nodes, A)) ;
113 :
114 66 : GRB_TRY (GrB_Matrix_new (&A_tril, type, num_nodes, num_nodes)) ;
115 66 : GRB_TRY (GrB_Matrix_new (&E, type, num_nodes, num_edges)) ;
116 66 : GRB_TRY (GrB_Matrix_new (&E_half, type, num_nodes, num_edges)) ;
117 :
118 : // get just the lower triangular entries
119 66 : GRB_TRY (GrB_select (A_tril, NULL, NULL, GrB_TRIL, A, 0, NULL)) ;
120 :
121 : #if GxB_IMPLEMENTATION < GxB_VERSION (10,0,0)
122 : bool is_uint64 = (type == GrB_UINT64) ;
123 : bool is_float = ((type == GrB_FP32) || (type == GrB_FP64)) ;
124 :
125 : // which intermediate type to cast to
126 : // uint64_t : 0, double: 1, int64_t: 2
127 : // if the input matrix type is GrB_INT* or GrB_UINT{8|16|32}, it becomes a int64_t
128 : // if the input matrix type is GrB_UINT64, it becomes a uint64_t
129 : // if the input matrix type is GrB_FP{32|64}, it becomes a double
130 : int which_type = (is_uint64 ? 0 : (is_float ? 1 : 2)) ;
131 :
132 : size_t value_sizes[3] = {sizeof(uint64_t), sizeof(double), sizeof(int64_t)} ;
133 :
134 : // (*result) = E ;
135 : // return (GrB_SUCCESS) ;
136 :
137 : // Arrays to extract A into
138 : LG_TRY (LAGraph_Malloc ((void**)(&row_indices), num_edges, sizeof(GrB_Index), msg)) ;
139 : LG_TRY (LAGraph_Malloc ((void**)(&col_indices), num_edges, sizeof(GrB_Index), msg)) ;
140 : LG_TRY (LAGraph_Malloc ((void**)(&values), num_edges, value_sizes[which_type], msg)) ;
141 :
142 : switch (which_type){
143 : case 0:
144 : GRB_TRY (GrB_Matrix_extractTuples_UINT64 (row_indices, col_indices, values, &num_edges, A_tril)) ;
145 : break;
146 : case 1:
147 : GRB_TRY (GrB_Matrix_extractTuples_FP64 (row_indices, col_indices, values, &num_edges, A_tril)) ;
148 : break;
149 : case 2:
150 : GRB_TRY (GrB_Matrix_extractTuples_INT64 (row_indices, col_indices, values, &num_edges, A_tril)) ;
151 : break;
152 : }
153 :
154 : #ifdef dbg
155 : printf("Printing A_tril values\n");
156 : for (int i = 0; i < num_edges; i++){
157 : printf("%ld %ld %.5f\n", row_indices[i], col_indices[i], values[i]);
158 : }
159 : #endif
160 : LG_TRY (LAGraph_Malloc ((void**)(&ramp), num_edges, sizeof(GrB_Index), msg)) ;
161 :
162 : #pragma omp parallel for
163 : for (GrB_Index i = 0 ; i < num_edges ; i++) {
164 : ramp[i] = i ;
165 : }
166 :
167 : // build E_1 with (row_indices, ramp, values)
168 : // build E_2 with (col_indices, ramp, values)
169 : // E = E_1 + E_2
170 : switch (which_type) {
171 : case 0:
172 : GRB_TRY (GrB_Matrix_build_UINT64 (E_half, col_indices, ramp, values, num_edges, NULL)) ;
173 : GRB_TRY (GrB_Matrix_build_UINT64 (E, row_indices, ramp, values, num_edges, NULL)) ;
174 : break;
175 : case 1:
176 : GRB_TRY (GrB_Matrix_build_FP64 (E_half, col_indices, ramp, values, num_edges, NULL)) ;
177 : GRB_TRY (GrB_Matrix_build_FP64 (E, row_indices, ramp, values, num_edges, NULL)) ;
178 : break;
179 : case 2:
180 : GRB_TRY (GrB_Matrix_build_INT64 (E_half, col_indices, ramp, values, num_edges, NULL)) ;
181 : GRB_TRY (GrB_Matrix_build_INT64 (E, row_indices, ramp, values, num_edges, NULL)) ;
182 : break;
183 : }
184 : GRB_TRY (GrB_eWiseAdd (E, NULL, NULL, GrB_PLUS_FP64, E, E_half, NULL)) ;
185 : #else
186 : // The vector types and sizes get overriden by extract.
187 66 : GRB_TRY (GrB_Vector_new(&x, GrB_BOOL, 0)) ;
188 66 : GRB_TRY (GrB_Vector_new(&i, GrB_BOOL, 0)) ;
189 66 : GRB_TRY (GrB_Vector_new(&j, GrB_BOOL, 0)) ;
190 66 : GRB_TRY (GxB_Matrix_extractTuples_Vector(i, j, x, A_tril, NULL)) ;
191 :
192 : // FUTURE: x is always size nvals(A), and never iso-valued, even if A_tril
193 : // is iso-valued. If A_tril is iso-valued, then Ex below could be
194 : // length 1 and con->iso could be false.
195 :
196 : // this load trick is quicker, but returns GrB_COLMAJOR
197 : #ifdef LOADTRICKIM
198 66 : GRB_TRY (GrB_Vector_new(&Ep, GrB_INT64, num_edges + 1)) ;
199 66 : GrB_Type ij_type = NULL;
200 66 : GRB_TRY (GxB_Vector_type(&ij_type, i));
201 66 : GRB_TRY (GrB_Vector_new(&Ei, ij_type, num_edges * 2)) ;
202 66 : GRB_TRY (GrB_assign(
203 : Ep, NULL, NULL, (int64_t) 1, GrB_ALL, 0, NULL));
204 :
205 : // Shuffle i and j into Ei.
206 : // Ei = [j[0], i[0], j[1], i[1], . . ., i[num_edges -1]]
207 : // Filling out Ei helps assign be much quicker.
208 66 : GRB_TRY (GrB_assign(
209 : Ei, NULL, NULL, (int64_t) 1, GrB_ALL, 0, NULL));
210 66 : GrB_Index stride[] = {(GrB_Index) 0, num_edges * 2 - 1, (GrB_Index) 2} ;
211 : // printf ("num_edges: %lu\n", num_edges) ;
212 : // printf ("stride: [%lu %lu %lu]\n", stride [0], stride [1], stride [2]) ;
213 66 : if (num_edges > 0)
214 : {
215 64 : GRB_TRY (GrB_Vector_assign(
216 : Ei, NULL, NULL, j, stride, GxB_STRIDE, NULL)) ;
217 64 : stride[GxB_BEGIN] = 1;
218 64 : GRB_TRY (GrB_Vector_assign(
219 : Ei, NULL, NULL, i, stride, GxB_STRIDE, NULL)) ;
220 : }
221 :
222 : // Ep = [0,2,4,...,2 * numedges]
223 66 : GRB_TRY (GrB_Vector_apply_IndexOp_INT64(
224 : Ep, NULL, NULL, GrB_ROWINDEX_INT64, Ep, (uint64_t) 0, NULL)) ;
225 66 : GRB_TRY (GrB_Vector_apply_BinaryOp2nd_INT64(
226 : Ep, NULL, NULL, GrB_TIMES_INT64, Ep, (int64_t) 2, NULL)) ;
227 :
228 : // Filling out Ex helps assign be much quicker.
229 66 : GRB_TRY (GrB_Vector_new(&Ex, type, num_edges * 2)) ;
230 66 : GRB_TRY (GrB_assign(
231 : Ex, NULL, NULL, (int64_t) 1, GrB_ALL, 0, NULL));
232 66 : stride[GxB_BEGIN] = 0;
233 66 : if (num_edges > 0)
234 : {
235 64 : GRB_TRY (GrB_Vector_assign(
236 : Ex, NULL, NULL, x, stride, GxB_STRIDE, NULL)) ;
237 64 : stride[GxB_BEGIN] = 1;
238 64 : GRB_TRY (GrB_Vector_assign(
239 : Ex, NULL, NULL, x, stride, GxB_STRIDE, NULL)) ;
240 : }
241 :
242 : //load up container
243 66 : GRB_TRY (GxB_Container_new(&con));
244 66 : GRB_TRY (GrB_free(&con->p));
245 66 : GRB_TRY (GrB_free(&con->i));
246 66 : GRB_TRY (GrB_free(&con->x));
247 66 : con->p = Ep; Ep = NULL ;
248 66 : con->i = Ei; Ei = NULL ;
249 66 : con->x = Ex; Ex = NULL ;
250 66 : con->format = GxB_SPARSE;
251 66 : con->orientation = GrB_COLMAJOR;
252 66 : con->nrows = num_nodes;
253 66 : con->ncols = num_edges;
254 66 : con->nvals = num_edges * 2;
255 66 : con->jumbled = false;
256 66 : con->iso = false;
257 : // Ep = [0,2,4,...,2 * numedges]
258 : // Ex = [x[0], x[0], x[1], x[1], . . ., x[num_edges -1]]
259 : // Ei = [j[0], i[0], j[1], i[1], . . ., i[num_edges -1]]
260 : // So each column k has two entries at j[k] and i[k] with values x[k]
261 66 : GRB_TRY (GxB_load_Matrix_from_Container(E, con, NULL));
262 66 : GRB_TRY (GrB_free(&con));
263 : #else
264 : GRB_TRY (GrB_Vector_new(
265 : &fullx, GrB_BOOL, num_edges)) ;
266 : GRB_TRY (GrB_assign(
267 : fullx, NULL, NULL, (bool) 1, GrB_ALL, 0, NULL));
268 : GRB_TRY (GrB_Descriptor_new(&build_desc));
269 : GRB_TRY (GrB_set(build_desc, GxB_USE_INDICES, GxB_COLINDEX_LIST));
270 : // fullx interpreted by index so is just a ramp.
271 : GRB_TRY (GxB_Matrix_build_Vector(E_half, j, fullx, x, NULL, build_desc));
272 : GRB_TRY (GxB_Matrix_build_Vector(E, i, fullx, x, NULL, build_desc));
273 : GRB_TRY (GrB_eWiseAdd (E, NULL, NULL, GrB_PLUS_FP64, E, E_half, NULL)) ;
274 : #endif
275 : #endif
276 : // LAGraph_Matrix_Print (E, LAGraph_COMPLETE, stdout, msg) ;
277 66 : LG_FREE_ALL ;
278 :
279 66 : (*result) = E ;
280 66 : return (GrB_SUCCESS) ;
281 : #else
282 : return (GrB_NOT_IMPLEMENTED) ;
283 : #endif
284 : }
|
