Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_msf.c
3 : //------------------------------------------------------------------------------
4 :
5 : // LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved.
6 : // SPDX-License-Identifier: BSD-2-Clause
7 : //
8 : // For additional details (including references to third party source code and
9 : // other files) see the LICENSE file or contact permission@sei.cmu.edu. See
10 : // Contributors.txt for a full list of contributors. Created, in part, with
11 : // funding and support from the U.S. Government (see Acknowledgments.txt file).
12 : // DM22-0790
13 :
14 : // Contributed by Yongzhe Zhang (zyz915@gmail.com)
15 : // Revised by Gabriel Gomez and Tim Davis
16 :
17 : //------------------------------------------------------------------------------
18 :
19 : /**
20 : * Code is based on Boruvka's minimum spanning forest algorithm
21 : *
22 : * The algorithm calculates the minimum spanning forest of a graph represented
23 : * by matrix A. It outputs a matrix containing lowest weight edges which connect
24 : * the nodes in the forest, and a vector indicating the connected component
25 : * "tree" which each node belongs to.
26 : */
27 :
28 : // TODO: not ready for src but getting close.
29 :
30 : // TODO: a "sanitize" input is fine for now in the experimental folder, but it
31 : // doesn't fit with the standard LAGraph API. It will need to be removed when
32 : // this method is moved to the src folder. The input will also become an
33 : // LAGraph_Graph, not a plain GrB_Matrix A.
34 :
35 : // FUTURE: the idx type could be INT32 if the matrix has < 2 billion nodes,
36 : // and the weight types could be expanded (int32 and float).
37 :
38 : #include "LG_internal.h"
39 : #include <LAGraph.h>
40 : #include <LAGraphX.h>
41 :
42 : //------------------------------------------------------------------------------
43 : // tuple: a tuple containing (weight,index)
44 : //------------------------------------------------------------------------------
45 :
46 : // LG_MSF_tuple_int is used if the input graph uses integer weights of any type;
47 : // LG_MSF_tuple_fp is used if the input graph is FP32 or FP64. Likewise for the
48 : // other *_int and *_fp types and operators.
49 :
50 : typedef struct
51 : {
52 : int64_t wInt;
53 : uint64_t idx;
54 : } LG_MSF_tuple_int;
55 :
56 : #define TUPLE_INT \
57 : "typedef struct \n" \
58 : "{ \n" \
59 : " int64_t wInt; \n" \
60 : " uint64_t idx; \n" \
61 : "} LG_MSF_tuple_int;"
62 :
63 : typedef struct
64 : {
65 : double wFp;
66 : uint64_t idx;
67 : } LG_MSF_tuple_fp;
68 :
69 : #define TUPLE_FP \
70 : "typedef struct \n" \
71 : "{ \n" \
72 : " double wFp; \n" \
73 : " uint64_t idx; \n" \
74 : "} LG_MSF_tuple_fp;"
75 :
76 : //------------------------------------------------------------------------------
77 : // context_type: context for IndexUnaryOps (using the theta input)
78 : //------------------------------------------------------------------------------
79 :
80 : typedef struct
81 : {
82 : uint64_t *parent; // parent of each vertex in the spanning forest
83 : struct
84 : {
85 : int64_t wInt;
86 : uint64_t idx;
87 : } *w_partner; // partner vertex in the spanning forest
88 : } LG_MSF_context_int;
89 :
90 : #define LG_MSF_CONTEXT_INT \
91 : "typedef struct \n" \
92 : "{ \n" \
93 : " uint64_t *parent; \n" \
94 : " struct \n" \
95 : " { \n" \
96 : " int64_t wInt; \n" \
97 : " uint64_t idx; \n" \
98 : " } *w_partner; \n" \
99 : "} LG_MSF_context_int;"
100 :
101 : typedef struct
102 : {
103 : uint64_t *parent; // parent of each vertex in the spanning forest
104 : struct
105 : {
106 : double wFp;
107 : uint64_t idx;
108 : } *w_partner; // partner vertex in the spanning forest
109 : } LG_MSF_context_fp;
110 :
111 : #define LG_MSF_CONTEXT_FP \
112 : "typedef struct \n" \
113 : "{ \n" \
114 : " uint64_t *parent; \n" \
115 : " struct \n" \
116 : " { \n" \
117 : " double wFp; \n" \
118 : " uint64_t idx; \n" \
119 : " } *w_partner; \n" \
120 : "} LG_MSF_context_fp;"
121 :
122 : //------------------------------------------------------------------------------
123 : // selectEdge: index-unary operator to select edges of min weight
124 : //------------------------------------------------------------------------------
125 :
126 : // generate solution:
127 : // for each element A(i, j), it is selected if
128 : // 1. weight[i] == A(i, j) -- where weight[i] stores i's minimum edge weight
129 : // 2. parent[j] == partner[i] -- j belongs to the specified connected component
130 :
131 16028 : void LG_MSF_selectEdge_int
132 : (
133 : bool *z,
134 : const int64_t *x,
135 : GrB_Index i,
136 : GrB_Index j,
137 : const LG_MSF_context_int *theta
138 : )
139 : {
140 32033 : (*z) = (theta->w_partner[i].wInt == *x) &&
141 16005 : (theta->parent[j] == theta->w_partner[i].idx);
142 16028 : }
143 :
144 : #define SELECTEDGE_INT \
145 : "void LG_MSF_selectEdge_int \n"\
146 : "( \n"\
147 : " bool *z, \n"\
148 : " const int64_t *x, \n"\
149 : " GrB_Index i, \n"\
150 : " GrB_Index j, \n"\
151 : " const LG_MSF_context_int *theta \n"\
152 : ") \n"\
153 : "{ \n"\
154 : " (*z) = (theta->w_partner[i].wInt == *x) && \n"\
155 : " (theta->parent[j] == theta->w_partner[i].idx);\n"\
156 : "}"
157 :
158 576770 : void LG_MSF_selectEdge_fp
159 : (
160 : bool *z,
161 : const double *x,
162 : GrB_Index i,
163 : GrB_Index j,
164 : const LG_MSF_context_fp *theta
165 : )
166 : {
167 598783 : (*z) = (theta->w_partner[i].wFp == *x) &&
168 22013 : (theta->parent[j] == theta->w_partner[i].idx);
169 576770 : }
170 :
171 : #define SELECTEDGE_FP \
172 : "void LG_MSF_selectEdge_fp \n"\
173 : "( \n"\
174 : " bool *z, \n"\
175 : " const double *x, \n"\
176 : " GrB_Index i, \n"\
177 : " GrB_Index j, \n"\
178 : " const LG_MSF_context_fp *theta \n"\
179 : ") \n"\
180 : "{ \n"\
181 : " (*z) = (theta->w_partner[i].wFp == *x) && \n"\
182 : " (theta->parent[j] == theta->w_partner[i].idx);\n"\
183 : "}"
184 :
185 : //------------------------------------------------------------------------------
186 : // removeEdge: remove edge (i,j) when i and j have the same parent
187 : //------------------------------------------------------------------------------
188 :
189 : // edge removal:
190 : // A(i, j) is removed when parent[i] == parent[j]
191 :
192 16028 : void LG_MSF_removeEdge_int
193 : (
194 : bool *z,
195 : const int64_t *x,
196 : GrB_Index i,
197 : GrB_Index j,
198 : const LG_MSF_context_int *theta
199 : )
200 : {
201 16028 : (*z) = (theta->parent[i] != theta->parent[j]);
202 16028 : }
203 :
204 : #define REMOVEEDGE_INT \
205 : "void LG_MSF_removeEdge_int \n"\
206 : "( \n"\
207 : " bool *z, \n"\
208 : " const int64_t *x, \n"\
209 : " GrB_Index i, \n"\
210 : " GrB_Index j, \n"\
211 : " const LG_MSF_context_int *theta \n"\
212 : ") \n"\
213 : "{ \n"\
214 : " (*z) = (theta->parent[i] != theta->parent[j]);\n"\
215 : "}"
216 :
217 576770 : void LG_MSF_removeEdge_fp
218 : (
219 : bool *z,
220 : const double *x,
221 : GrB_Index i,
222 : GrB_Index j,
223 : const LG_MSF_context_fp *theta
224 : )
225 : {
226 576770 : (*z) = (theta->parent[i] != theta->parent[j]);
227 576770 : }
228 :
229 : #define REMOVEEDGE_FP \
230 : "void LG_MSF_removeEdge_fp \n"\
231 : "( \n"\
232 : " bool *z, \n"\
233 : " const double *x, \n"\
234 : " GrB_Index i, \n"\
235 : " GrB_Index j, \n"\
236 : " const LG_MSF_context_fp *theta \n"\
237 : ") \n"\
238 : "{ \n"\
239 : " (*z) = (theta->parent[i] != theta->parent[j]);\n"\
240 : "}"
241 :
242 : //------------------------------------------------------------------------------
243 : // combine: create a tuple from a weight and an index
244 : //------------------------------------------------------------------------------
245 :
246 11750 : void LG_MSF_combine_int
247 : (
248 : LG_MSF_tuple_int *z,
249 : const int64_t *x,
250 : const uint64_t *y
251 : )
252 : {
253 11750 : z->wInt = *x;
254 11750 : z->idx = *y;
255 11750 : }
256 :
257 : #define COMBINE_INT \
258 : "void LG_MSF_combine_int \n"\
259 : "( \n"\
260 : " LG_MSF_tuple_int *z, \n"\
261 : " const int64_t *x, \n"\
262 : " const uint64_t *y \n"\
263 : ") \n"\
264 : "{ \n"\
265 : " z->wInt = *x; \n"\
266 : " z->idx = *y; \n"\
267 : "}"
268 :
269 300726 : void LG_MSF_combine_fp
270 : (
271 : LG_MSF_tuple_fp *z,
272 : const double *x,
273 : const uint64_t *y
274 : )
275 : {
276 300726 : z->wFp = *x;
277 300726 : z->idx = *y;
278 300726 : }
279 :
280 : #define COMBINE_FP \
281 : "void LG_MSF_combine_fp \n"\
282 : "( \n"\
283 : " LG_MSF_tuple_fp *z, \n"\
284 : " const double *x, \n"\
285 : " const uint64_t *y \n"\
286 : ") \n"\
287 : "{ \n"\
288 : " z->wFp = *x; \n"\
289 : " z->idx = *y; \n"\
290 : "}"
291 :
292 : //------------------------------------------------------------------------------
293 : // get_first: get first item in a tuple (the weight)
294 : //------------------------------------------------------------------------------
295 :
296 1220 : void LG_MSF_get_first_int (int64_t *y, const LG_MSF_tuple_int *x)
297 : {
298 1220 : *y = x->wInt;
299 1220 : }
300 :
301 : #define GET_FIRST_INT \
302 : "void LG_MSF_get_first_int (int64_t *y, const LG_MSF_tuple_int *x) \n" \
303 : "{ \n" \
304 : " *y = x->wInt; \n" \
305 : "}"
306 :
307 2076 : void LG_MSF_get_first_fp (double *y, const LG_MSF_tuple_fp *x)
308 : {
309 2076 : *y = x->wFp;
310 2076 : }
311 :
312 : #define GET_FIRST_FP \
313 : "void LG_MSF_get_first_fp (double *y, const LG_MSF_tuple_fp *x) \n" \
314 : "{ \n" \
315 : " *y = x->wFp; \n" \
316 : "}"
317 :
318 : //------------------------------------------------------------------------------
319 : // get_second: get second item in a tuple (the index)
320 : //------------------------------------------------------------------------------
321 :
322 3639 : void LG_MSF_get_second_int (uint64_t *y, const LG_MSF_tuple_int *x)
323 : {
324 3639 : *y = x->idx;
325 3639 : }
326 :
327 : #define GET_SECOND_INT \
328 : "void LG_MSF_get_second_int (uint64_t *y, const LG_MSF_tuple_int *x) \n" \
329 : "{ \n" \
330 : " *y = x->idx; \n" \
331 : "}"
332 :
333 12310 : void LG_MSF_get_second_fp (uint64_t *y, const LG_MSF_tuple_fp *x)
334 : {
335 12310 : *y = x->idx;
336 12310 : }
337 :
338 : #define GET_SECOND_FP \
339 : "void LG_MSF_get_second_fp (uint64_t *y, const LG_MSF_tuple_fp *x) \n" \
340 : "{ \n" \
341 : " *y = x->idx; \n" \
342 : "}"
343 :
344 : //------------------------------------------------------------------------------
345 : // tupleMin: z = the min tuple of x and y
346 : //------------------------------------------------------------------------------
347 :
348 10530 : void LG_MSF_tupleMin_int
349 : (
350 : LG_MSF_tuple_int *z,
351 : const LG_MSF_tuple_int *x,
352 : const LG_MSF_tuple_int *y
353 : )
354 : {
355 20613 : bool xSmaller = x->wInt < y->wInt ||
356 10083 : (x->wInt == y->wInt && x->idx < y->idx);
357 10530 : z->wInt = (xSmaller)? x->wInt: y->wInt;
358 10530 : z->idx = (xSmaller)? x->idx: y->idx;
359 10530 : }
360 :
361 : #define TUPLEMIN_INT \
362 : "void LG_MSF_tupleMin_int \n"\
363 : "( \n"\
364 : " LG_MSF_tuple_int *z, \n"\
365 : " const LG_MSF_tuple_int *x, \n"\
366 : " const LG_MSF_tuple_int *y \n"\
367 : ") \n"\
368 : "{ \n"\
369 : " bool xSmaller = x->wInt < y->wInt || \n"\
370 : " (x->wInt == y->wInt && x->idx < y->idx);\n"\
371 : " z->wInt = (xSmaller)? x->wInt: y->wInt; \n"\
372 : " z->idx = (xSmaller)? x->idx: y->idx; \n"\
373 : "}"
374 :
375 298650 : void LG_MSF_tupleMin_fp
376 : (
377 : LG_MSF_tuple_fp *z,
378 : const LG_MSF_tuple_fp *x,
379 : const LG_MSF_tuple_fp *y
380 : )
381 : {
382 298650 : bool xSmaller = x->wFp < y->wFp || (x->wFp == y->wFp && x->idx < y->idx);
383 298650 : z->wFp = (xSmaller)? x->wFp: y->wFp;
384 298650 : z->idx = (xSmaller)? x->idx: y->idx;
385 298650 : }
386 :
387 : #define TUPLEMIN_FP \
388 : "void LG_MSF_tupleMin_fp \n"\
389 : "( \n"\
390 : " LG_MSF_tuple_fp *z, \n"\
391 : " const LG_MSF_tuple_fp *x, \n"\
392 : " const LG_MSF_tuple_fp *y \n"\
393 : ") \n"\
394 : "{ \n"\
395 : " bool xSmaller = x->wFp < y->wFp || (x->wFp == y->wFp && x->idx < y->idx);\n"\
396 : " z->wFp = (xSmaller)? x->wFp: y->wFp; \n"\
397 : " z->idx = (xSmaller)? x->idx: y->idx; \n"\
398 : "}"
399 :
400 : //------------------------------------------------------------------------------
401 : // tuple2nd: z = y
402 : //------------------------------------------------------------------------------
403 :
404 2419 : void LG_MSF_tuple2nd_int
405 : (
406 : LG_MSF_tuple_int *z,
407 : const void *x,
408 : const LG_MSF_tuple_int *y
409 : )
410 : {
411 2419 : z->wInt = y->wInt;
412 2419 : z->idx = y->idx;
413 2419 : }
414 :
415 : #define TUPLE2ND_INT \
416 : "void LG_MSF_tuple2nd_int \n"\
417 : "( \n"\
418 : " LG_MSF_tuple_int *z, \n"\
419 : " const void *x, \n"\
420 : " const LG_MSF_tuple_int *y\n"\
421 : ") \n"\
422 : "{ \n"\
423 : " z->wInt = y->wInt; \n"\
424 : " z->idx = y->idx; \n"\
425 : "}"
426 :
427 10234 : void LG_MSF_tuple2nd_fp
428 : (
429 : LG_MSF_tuple_fp *z,
430 : const void *x,
431 : const LG_MSF_tuple_fp *y
432 : )
433 : {
434 10234 : z->wFp = y->wFp;
435 10234 : z->idx = y->idx;
436 10234 : }
437 :
438 : #define TUPLE2ND_FP \
439 : "void LG_MSF_tuple2nd_fp \n"\
440 : "( \n"\
441 : " LG_MSF_tuple_fp *z, \n"\
442 : " const void *x, \n"\
443 : " const LG_MSF_tuple_fp *y\n"\
444 : ") \n"\
445 : "{ \n"\
446 : " z->wFp = y->wFp; \n"\
447 : " z->idx = y->idx; \n"\
448 : "}"
449 :
450 : //------------------------------------------------------------------------------
451 : // tupleEq: true if two tuples are equal
452 : //------------------------------------------------------------------------------
453 :
454 2419 : void LG_MSF_tupleEq_int
455 : (
456 : bool *z,
457 : const LG_MSF_tuple_int *x,
458 : const LG_MSF_tuple_int *y
459 : )
460 : {
461 2419 : *z = (x->wInt == y->wInt) && (x->idx == y->idx);
462 2419 : }
463 :
464 : #define TUPLEEQ_INT \
465 : "void LG_MSF_tupleEq_int \n"\
466 : "( \n"\
467 : " bool *z, \n"\
468 : " const LG_MSF_tuple_int *x, \n"\
469 : " const LG_MSF_tuple_int *y \n"\
470 : ") \n"\
471 : "{ \n"\
472 : " *z = (x->wInt == y->wInt) && (x->idx == y->idx);\n"\
473 : "}"
474 :
475 10234 : void LG_MSF_tupleEq_fp
476 : (
477 : bool *z,
478 : const LG_MSF_tuple_fp *x,
479 : const LG_MSF_tuple_fp *y
480 : )
481 : {
482 10234 : *z = (x->wFp == y->wFp) && (x->idx == y->idx);
483 10234 : }
484 :
485 : #define TUPLEEQ_FP \
486 : "void LG_MSF_tupleEq_fp \n"\
487 : "( \n"\
488 : " bool *z, \n"\
489 : " const LG_MSF_tuple_fp *x, \n"\
490 : " const LG_MSF_tuple_fp *y \n"\
491 : ") \n"\
492 : "{ \n"\
493 : " *z = (x->wFp == y->wFp) && (x->idx == y->idx);\n"\
494 : "}"
495 :
496 : //------------------------------------------------------------------------------
497 :
498 : #undef LG_FREE_ALL
499 : #define LG_FREE_ALL \
500 : { \
501 : GrB_free (&S); \
502 : GrB_free (&T); \
503 : LAGraph_Free ((void **) &SI, msg); \
504 : LAGraph_Free ((void **) &SJ, msg); \
505 : LAGraph_Free ((void **) &SX, msg); \
506 : LAGraph_Free ((void **) &context_int.parent, msg); \
507 : LAGraph_Free ((void **) &context_fp.parent, msg); \
508 : GrB_free (&f); \
509 : GrB_free (&w); \
510 : GrB_free (&I); \
511 : GrB_free (&t); \
512 : GrB_free (&edge); \
513 : GrB_free (&cedge); \
514 : GrB_free (&tedge); \
515 : GrB_free (&mask); \
516 : GrB_free (&index_v); \
517 : GrB_free (&combine); \
518 : GrB_free (&minComb); \
519 : GrB_free (&get_first); \
520 : GrB_free (&get_second); \
521 : GrB_free (&selectEdge); \
522 : GrB_free (&removeEdge); \
523 : GrB_free (&context_type); \
524 : GrB_free (&parent_v); \
525 : GrB_free (&ramp); \
526 : GrB_free (&tupleMin); \
527 : GrB_free (&tuple2nd); \
528 : GrB_free (&tupleEq); \
529 : GrB_free (&tupleMin_monoid); \
530 : GrB_free (&tupleMin2nd); \
531 : GrB_free (&tuple); \
532 : GrB_free (&max_weight); \
533 : }
534 :
535 : #ifdef DBG
536 : #undef GRB_CATCH
537 : #define GRB_CATCH(info) \
538 : { \
539 : printf ("GraphBLAS failure (file %s, line %d): info: %d", \
540 : __FILE__, __LINE__, info) ; \
541 : LG_ERROR_MSG ("GraphBLAS failure (file %s, line %d): info: %d", \
542 : __FILE__, __LINE__, info) ; \
543 : LG_FREE_ALL ; \
544 : return (info) ; \
545 : }
546 : #endif
547 :
548 : //------------------------------------------------------------------------------
549 : // dump_tuple_vector: debugging only
550 : //------------------------------------------------------------------------------
551 :
552 : // #define DBG
553 :
554 : #ifdef DBG
555 : static void dump_tuple_vector
556 : (
557 : char *vname,
558 : GrB_Vector v, // of type tuple (int or fp)
559 : GrB_Type weight_type
560 : )
561 : {
562 : GrB_Index n = 0 ;
563 : GrB_Info info ;
564 : GrB_Vector_size (&n, v) ;
565 : printf ("\ntuple vector %s, size %lu\n", vname, n) ;
566 : if (weight_type == GrB_INT64)
567 : {
568 : printf ("weight type: int64\n") ;
569 : LG_MSF_tuple_int e ;
570 : for (int i = 0 ; i < n ; i++)
571 : {
572 : info = GrB_Vector_extractElement_UDT (&e, v, i) ;
573 : if (info == GrB_SUCCESS)
574 : {
575 : printf (" (%d) (%ld, %lu)\n", i, e.wInt, e.idx) ;
576 : }
577 : }
578 : }
579 : else
580 : {
581 : printf ("weight type: double\n") ;
582 : LG_MSF_tuple_fp e ;
583 : for (int i = 0 ; i < n ; i++)
584 : {
585 : info = GrB_Vector_extractElement_UDT (&e, v, i) ;
586 : if (info == GrB_SUCCESS)
587 : {
588 : printf (" (%d) (%g, %lu)\n", i, e.wFp, e.idx) ;
589 : }
590 : }
591 : }
592 : }
593 : #endif
594 :
595 : //------------------------------------------------------------------------------
596 : // LAGraph_msf
597 : //------------------------------------------------------------------------------
598 :
599 29 : int LAGraph_msf
600 : (
601 : GrB_Matrix *forest_edges, // output: an unsymmetrical matrix, containing
602 : // the edges in the spanning forest
603 : GrB_Vector *componentId, // output: The connected component of each node
604 : // componentId[i] is the representative node of the
605 : // component set that i is in.
606 : GrB_Matrix A, // input matrix
607 : bool sanitize, // if true, ensure A is symmetric
608 : char *msg
609 : )
610 : {
611 : #if LG_SUITESPARSE_GRAPHBLAS_V10
612 29 : LG_CLEAR_MSG ;
613 :
614 29 : LG_MSF_context_int context_int = {.parent = NULL, .w_partner = NULL } ;
615 29 : LG_MSF_context_fp context_fp = {.parent = NULL, .w_partner = NULL } ;
616 29 : LG_MSF_tuple_int inf_int = {.wInt = INT64_MAX, .idx = UINT64_MAX};
617 29 : LG_MSF_tuple_fp inf_fp = {.wFp = INFINITY , .idx = UINT64_MAX};
618 :
619 : GrB_Info info;
620 : GrB_Index n;
621 29 : GrB_Matrix S = NULL, T = NULL;
622 29 : GrB_Vector f = NULL, I = NULL, t = NULL, parent_v = NULL, tedge = NULL,
623 29 : edge = NULL, cedge = NULL, mask = NULL, index_v = NULL, ramp = NULL,
624 29 : w = NULL ;
625 :
626 29 : GrB_Index *SI = NULL, *SJ = NULL;
627 29 : void *SX = NULL;
628 29 : GrB_Type context_type = NULL, tuple = NULL, weight_type = NULL,
629 29 : ignore = NULL ;
630 29 : GrB_BinaryOp combine = NULL, tupleMin = NULL, tuple2nd = NULL,
631 29 : tupleEq = NULL;
632 29 : GrB_Monoid tupleMin_monoid = NULL;
633 29 : GrB_Semiring minComb = NULL, tupleMin2nd = NULL;
634 29 : GrB_UnaryOp get_first = NULL, get_second = NULL;
635 29 : GrB_Scalar max_weight = NULL;
636 29 : int edge_handling = GrB_DEFAULT;
637 29 : uint64_t edge_size = 0, edge_n = 0;
638 29 : GrB_IndexUnaryOp selectEdge = NULL, removeEdge = NULL;
639 :
640 : //--------------------------------------------------------------------------
641 : // Check inputs
642 : //--------------------------------------------------------------------------
643 :
644 29 : if (forest_edges == NULL || A == NULL) return (GrB_NULL_POINTER) ;
645 : GrB_Index ncols ;
646 28 : GRB_TRY (GrB_Matrix_nrows (&n, A)) ;
647 28 : GRB_TRY (GrB_Matrix_ncols (&ncols, A)) ;
648 28 : LG_ASSERT(n == ncols, GrB_DIMENSION_MISMATCH) ;
649 :
650 27 : GrB_BinaryOp min_weight = NULL ;
651 27 : size_t sx_size = 0 ;
652 :
653 27 : GrB_Type_Code tcode = 0;
654 27 : GrB_Matrix_get_INT32(A, (int *) &(tcode), GrB_EL_TYPE_CODE);
655 27 : switch (tcode)
656 : {
657 19 : case GrB_INT8_CODE:
658 : case GrB_INT16_CODE:
659 : case GrB_INT32_CODE:
660 : case GrB_INT64_CODE:
661 : case GrB_BOOL_CODE:
662 : case GrB_UINT8_CODE:
663 : case GrB_UINT16_CODE:
664 : case GrB_UINT32_CODE:
665 : case GrB_UINT64_CODE:
666 : // integer edge weights: use INT64
667 19 : weight_type = GrB_INT64 ;
668 19 : min_weight = GrB_MIN_INT64 ;
669 19 : sx_size = sizeof (int64_t) ;
670 19 : break;
671 :
672 7 : case GrB_FP32_CODE:
673 : case GrB_FP64_CODE:
674 : // floating-point edge weights: use FP64
675 7 : weight_type = GrB_FP64 ;
676 7 : min_weight = GrB_MIN_FP64 ;
677 7 : sx_size = sizeof (double) ;
678 7 : break;
679 :
680 1 : default:
681 : // other types are not supported
682 1 : LG_ASSERT(false, GrB_DOMAIN_MISMATCH) ;
683 : break;
684 : }
685 :
686 : //--------------------------------------------------------------------------
687 : // create types and operators
688 : //--------------------------------------------------------------------------
689 :
690 26 : GRB_TRY (GxB_Scalar_new(&max_weight, weight_type)) ;
691 26 : void *inf = NULL ;
692 :
693 26 : if (weight_type == GrB_INT64)
694 : {
695 :
696 : //----------------------------------------------------------------------
697 : // types and ops for INT64 weights
698 : //----------------------------------------------------------------------
699 :
700 19 : GRB_TRY (GxB_Type_new (&tuple, sizeof (LG_MSF_tuple_int),
701 : "LG_MSF_tuple_int", TUPLE_INT)) ;
702 :
703 19 : GRB_TRY (GxB_BinaryOp_new (
704 : &combine, (GxB_binary_function) LG_MSF_combine_int,
705 : tuple, weight_type, GrB_UINT64,
706 : "LG_MSF_combine_int", COMBINE_INT)) ;
707 :
708 19 : GRB_TRY (GxB_Scalar_setElement_INT64(max_weight, INT64_MAX)) ;
709 :
710 19 : GRB_TRY (GxB_BinaryOp_new (
711 : &tupleMin, (GxB_binary_function) LG_MSF_tupleMin_int,
712 : tuple, tuple, tuple,
713 : "LG_MSF_tupleMin_int", TUPLEMIN_INT)) ;
714 :
715 19 : GRB_TRY (GxB_BinaryOp_new (
716 : &tuple2nd, (GxB_binary_function) LG_MSF_tuple2nd_int,
717 : tuple, GrB_BOOL, tuple,
718 : "LG_MSF_tuple2nd_int", TUPLE2ND_INT)) ;
719 :
720 19 : GRB_TRY (GxB_BinaryOp_new (
721 : &tupleEq, (GxB_binary_function) LG_MSF_tupleEq_int,
722 : GrB_BOOL, tuple, tuple,
723 : "LG_MSF_tupleEq_int", TUPLEEQ_INT)) ;
724 :
725 19 : inf = (void *) (&inf_int) ;
726 :
727 19 : GRB_TRY (GxB_UnaryOp_new (
728 : &get_first, (GxB_unary_function) LG_MSF_get_first_int, weight_type,
729 : tuple, "LG_MSF_get_first_int", GET_FIRST_INT)) ;
730 :
731 19 : GRB_TRY (GxB_UnaryOp_new (
732 : &get_second, (GxB_unary_function) LG_MSF_get_second_int, GrB_UINT64,
733 : tuple, "LG_MSF_get_second_int", GET_SECOND_INT)) ;
734 :
735 : // context type
736 19 : GRB_TRY (GxB_Type_new (
737 : &context_type, sizeof (LG_MSF_context_int),
738 : "LG_MSF_context_int", LG_MSF_CONTEXT_INT)) ;
739 :
740 : // ops for GrB_select
741 19 : GRB_TRY(GxB_IndexUnaryOp_new (
742 : &selectEdge, (GxB_index_unary_function) LG_MSF_selectEdge_int,
743 : GrB_BOOL, weight_type, context_type,
744 : "LG_MSF_selectEdge_int", SELECTEDGE_INT)) ;
745 :
746 19 : GRB_TRY(GxB_IndexUnaryOp_new (
747 : &removeEdge, (void *) LG_MSF_removeEdge_int, GrB_BOOL, weight_type,
748 : context_type, "LG_MSF_removeEdge_int", REMOVEEDGE_INT)) ;
749 :
750 : }
751 : else
752 : {
753 :
754 : //----------------------------------------------------------------------
755 : // types and ops for FP64 weights
756 : //----------------------------------------------------------------------
757 :
758 7 : GRB_TRY (GxB_Type_new (&tuple, sizeof (LG_MSF_tuple_fp),
759 : "LG_MSF_tuple_fp", TUPLE_FP)) ;
760 :
761 7 : GRB_TRY (GxB_BinaryOp_new (
762 : &combine, (GxB_binary_function) LG_MSF_combine_fp,
763 : tuple, weight_type, GrB_UINT64,
764 : "LG_MSF_combine_fp", COMBINE_FP)) ;
765 :
766 7 : GRB_TRY (GxB_Scalar_setElement_FP64(max_weight, INFINITY)) ;
767 :
768 7 : GRB_TRY (GxB_BinaryOp_new (
769 : &tupleMin, (GxB_binary_function) LG_MSF_tupleMin_fp,
770 : tuple, tuple, tuple,
771 : "LG_MSF_tupleMin_fp", TUPLEMIN_FP)) ;
772 :
773 7 : GRB_TRY (GxB_BinaryOp_new (
774 : &tuple2nd, (GxB_binary_function) LG_MSF_tuple2nd_fp,
775 : tuple, GrB_BOOL, tuple,
776 : "LG_MSF_tuple2nd_fp", TUPLE2ND_FP)) ;
777 :
778 7 : GRB_TRY (GxB_BinaryOp_new (
779 : &tupleEq, (GxB_binary_function) LG_MSF_tupleEq_fp,
780 : GrB_BOOL, tuple, tuple,
781 : "LG_MSF_tupleEq_fp", TUPLEEQ_FP)) ;
782 :
783 7 : inf = (void *) (&inf_fp) ;
784 :
785 7 : GRB_TRY (GxB_UnaryOp_new (
786 : &get_first, (GxB_unary_function) LG_MSF_get_first_fp, weight_type,
787 : tuple, "LG_MSF_get_first_fp", GET_FIRST_FP)) ;
788 :
789 7 : GRB_TRY (GxB_UnaryOp_new (
790 : &get_second, (GxB_unary_function) LG_MSF_get_second_fp, GrB_UINT64,
791 : tuple, "LG_MSF_get_second_fp", GET_SECOND_FP)) ;
792 :
793 7 : GRB_TRY (GxB_Type_new (
794 : &context_type, sizeof (LG_MSF_context_fp),
795 : "LG_MSF_context_fp", LG_MSF_CONTEXT_FP)) ;
796 :
797 : // ops for GrB_select
798 7 : GRB_TRY(GxB_IndexUnaryOp_new (
799 : &selectEdge, (GxB_index_unary_function) LG_MSF_selectEdge_fp,
800 : GrB_BOOL, weight_type, context_type,
801 : "LG_MSF_selectEdge_fp", SELECTEDGE_FP)) ;
802 :
803 7 : GRB_TRY(GxB_IndexUnaryOp_new (
804 : &removeEdge, (void *) LG_MSF_removeEdge_fp, GrB_BOOL, weight_type,
805 : context_type, "LG_MSF_removeEdge_fp", REMOVEEDGE_FP)) ;
806 : }
807 :
808 26 : GRB_TRY (GrB_Monoid_new_UDT (&tupleMin_monoid, tupleMin, inf)) ;
809 26 : GRB_TRY (GrB_Semiring_new (&minComb, tupleMin_monoid, combine)) ;
810 26 : GRB_TRY (GrB_Semiring_new (&tupleMin2nd, tupleMin_monoid, tuple2nd)) ;
811 :
812 : //--------------------------------------------------------------------------
813 : // create matrices and vectors
814 : //--------------------------------------------------------------------------
815 :
816 26 : GRB_TRY (GrB_Matrix_new (&S, weight_type, n, n)) ;
817 :
818 26 : if (sanitize)
819 : {
820 : // S = A+A', and typecasting to weight_type
821 6 : GRB_TRY (GrB_Matrix_eWiseAdd_BinaryOp
822 : (S, NULL, NULL, min_weight, A, A, GrB_DESC_T1)) ;
823 : }
824 : else
825 : {
826 : // S = A, typecasting to weight_type, if necessary
827 20 : GRB_TRY (GrB_Matrix_assign
828 : (S, NULL, NULL, A, GrB_ALL, n, GrB_ALL, n, NULL)) ;
829 : }
830 :
831 26 : GRB_TRY (GrB_Matrix_new (&T, weight_type, n, n)) ;
832 26 : GRB_TRY (GrB_Vector_new (&w, weight_type, n)) ;
833 26 : GRB_TRY (GrB_Vector_new (&t, GrB_UINT64, n)) ;
834 26 : GRB_TRY (GrB_Vector_new (&f, GrB_UINT64, n)) ;
835 26 : GRB_TRY (GrB_Vector_new (&ramp, GrB_INT64, n + 1)) ;
836 26 : GRB_TRY (GrB_Vector_new (&edge, tuple, n)) ;
837 26 : GRB_TRY (GrB_Vector_new (&cedge, tuple, n)) ;
838 26 : GRB_TRY (GrB_Vector_new (&tedge, tuple, n)) ;
839 26 : GRB_TRY (GrB_Vector_new (&mask, GrB_BOOL, n)) ;
840 26 : GRB_TRY (GrB_Vector_new (&index_v, GrB_UINT64, n)) ;
841 26 : GRB_TRY (GrB_Vector_new (&parent_v, GrB_UINT64, n)) ;
842 :
843 26 : LG_TRY (LAGraph_Malloc ((void **) &SI, 2*n, sizeof (GrB_Index), msg)) ;
844 26 : LG_TRY (LAGraph_Malloc ((void **) &SJ, 2*n, sizeof (GrB_Index), msg)) ;
845 26 : LG_TRY (LAGraph_Malloc (&SX, 2*n, sx_size, msg)) ;
846 :
847 : // prepare vectors
848 26 : GRB_TRY (GrB_Vector_assign_UINT64 (
849 : f, NULL, NULL, (uint64_t) 0, GrB_ALL, n, NULL)) ;
850 26 : GRB_TRY (GrB_Vector_apply_IndexOp_INT64 (
851 : f, NULL, NULL, GrB_ROWINDEX_INT64, f, (int64_t) 0, NULL)) ;
852 26 : GRB_TRY (GrB_Vector_dup (&I, f)) ;
853 26 : GRB_TRY (GrB_Vector_assign_UINT64 (
854 : ramp, NULL, NULL, (uint64_t) 0, GrB_ALL, n + 1, NULL)) ;
855 26 : GRB_TRY (GrB_Vector_apply_IndexOp_INT64 (
856 : ramp, NULL, NULL, GrB_ROWINDEX_INT64, ramp, (int64_t) 0, NULL)) ;
857 :
858 : // create context
859 26 : if (weight_type == GrB_INT64)
860 : {
861 19 : LG_TRY (LAGraph_Malloc
862 : ((void **) &context_int.parent, n, sizeof (uint64_t), msg)) ;
863 19 : GRB_TRY (GrB_Vector_extractTuples (NULL, context_int.parent, &n, f)) ;
864 19 : GRB_TRY (GxB_Vector_load(parent_v, (void **) &context_int.parent,
865 : GrB_UINT64, n, n * sizeof (uint64_t), GxB_IS_READONLY, NULL)) ;
866 : }
867 : else
868 : {
869 7 : LG_TRY (LAGraph_Malloc
870 : ((void **) &context_fp.parent, n, sizeof (double), msg)) ;
871 7 : GRB_TRY (GrB_Vector_extractTuples (NULL, context_fp.parent, &n, f)) ;
872 7 : GRB_TRY (GxB_Vector_load(parent_v, (void **) &context_fp.parent,
873 : GrB_UINT64, n, n * sizeof (double), GxB_IS_READONLY, NULL)) ;
874 : }
875 :
876 : //--------------------------------------------------------------------------
877 : // the main computation
878 : //--------------------------------------------------------------------------
879 :
880 26 : GrB_Index nvals, ntuples = 0, num;
881 26 : bool diff = false;
882 26 : GRB_TRY (GrB_Matrix_nvals (&nvals, S)) ;
883 76 : for (int iters = 1; nvals > 0; iters++)
884 : {
885 : // every vertex points to a root vertex at the beginning
886 : // edge[u] = u's minimum edge (weight and index are encoded together)
887 50 : GRB_TRY (GrB_Vector_assign_UDT (
888 : edge, NULL, NULL, inf, GrB_ALL, 0, NULL)) ;
889 :
890 : // each edge looks at its adjacent edges and picks the one with the
891 : // minimum weight. This weight is put into a tuple with the
892 : // representative value of the connected componect the edge connects to
893 50 : GRB_TRY (GrB_mxv (edge, NULL, tupleMin, minComb, S, f, NULL)) ;
894 :
895 : // cedge[u] = children's minimum edge | if u is a root
896 : // = (max_weight, u) | otherwise
897 50 : GRB_TRY (GrB_Vector_apply_BinaryOp1st_Scalar (
898 : cedge, NULL, NULL, combine, max_weight, I, NULL)) ;
899 50 : LG_TRY (LAGraph_FastAssign_Semiring(
900 : cedge, NULL, tupleMin, parent_v, edge, ramp, tupleMin2nd, NULL, msg
901 : )) ;
902 :
903 : // if (f[u] == u) f[u] := get_second(cedge[u]) -- the index
904 50 : GRB_TRY (GrB_eWiseMult (mask, NULL, NULL, GrB_EQ_UINT64, f, I, NULL)) ;
905 50 : GRB_TRY (GrB_apply (
906 : f, mask, GrB_SECOND_UINT64, get_second, cedge, NULL)) ;
907 :
908 : // identify all the vertex pairs (u, v) where f[u] == v and f[v] == u
909 : // and then select the minimum of u, v as the new root;
910 : // if (f[f[i]] == i) f[i] = min(f[i], i)
911 50 : GRB_TRY (GxB_Vector_extract_Vector (t, NULL, NULL, f, f, NULL)) ;
912 50 : GRB_TRY (GrB_eWiseMult (mask, NULL, NULL, GrB_EQ_UINT64, I, t, NULL)) ;
913 50 : GRB_TRY (GrB_assign (f, mask, GrB_MIN_UINT64, I, GrB_ALL, 0, NULL)) ;
914 :
915 : // five steps to generate the solution
916 : // 1. new roots (f[i] == i) revise their entries in cedge
917 50 : GRB_TRY (GrB_eWiseMult (mask, NULL, NULL, GrB_EQ_UINT64, I, f, NULL)) ;
918 50 : GRB_TRY (GrB_assign (cedge, mask, NULL, inf, GrB_ALL, 0, NULL)) ;
919 :
920 : // 2. every vertex tries to know whether one of its edges is selected
921 50 : GRB_TRY (GxB_Vector_extract_Vector (
922 : tedge, NULL, NULL, cedge, parent_v, NULL)) ;
923 50 : GRB_TRY (GrB_eWiseMult (mask ,NULL, NULL, tupleEq, edge, tedge, NULL)) ;
924 :
925 : // 3. each root picks a vertex from its children to
926 : // generate the solution
927 50 : GRB_TRY (GrB_assign (index_v, NULL, NULL, n, GrB_ALL, 0, NULL)) ;
928 50 : GRB_TRY (GrB_assign (index_v, mask, NULL, I, GrB_ALL, 0, NULL)) ;
929 50 : GRB_TRY (GrB_assign (t, NULL, NULL, n, GrB_ALL, 0, NULL)) ;
930 50 : LG_TRY (LAGraph_FastAssign_Semiring(
931 : t, NULL, GrB_MIN_UINT64, parent_v, index_v, ramp,
932 : GrB_MIN_SECOND_SEMIRING_UINT64, NULL, msg
933 : )) ;
934 50 : GRB_TRY (GxB_Vector_extract_Vector (
935 : index_v, NULL, NULL, t, parent_v, NULL)) ;
936 50 : GRB_TRY (GrB_eWiseMult (
937 : mask ,NULL, NULL, GrB_EQ_UINT64, I, index_v, NULL)) ;
938 :
939 : // 4. generate the select function
940 50 : if (weight_type == GrB_INT64)
941 : {
942 28 : GRB_TRY (GxB_Vector_unload(
943 : edge, (void **) &context_int.w_partner, &ignore, &edge_n,
944 : &edge_size, &edge_handling, NULL)) ;
945 28 : GRB_TRY (GrB_Matrix_select_UDT (T, NULL, NULL, selectEdge, S,
946 : &context_int, NULL)) ;
947 28 : GRB_TRY (GxB_Vector_load(
948 : edge, (void **) &context_int.w_partner, tuple, edge_n,
949 : edge_size, edge_handling, NULL)) ;
950 : }
951 : else
952 : {
953 22 : GRB_TRY (GxB_Vector_unload(
954 : edge, (void **) &context_fp.w_partner, &ignore, &edge_n,
955 : &edge_size, &edge_handling, NULL)) ;
956 22 : GRB_TRY (GrB_Matrix_select_UDT (
957 : T, NULL, NULL, selectEdge, S, &context_fp, NULL)) ;
958 22 : GRB_TRY (GxB_Vector_load(
959 : edge, (void **) &context_fp.w_partner, tuple, edge_n, edge_size,
960 : edge_handling, NULL)) ;
961 : }
962 :
963 50 : GRB_TRY (GrB_Vector_clear (t)) ;
964 :
965 : // 5. the generated matrix may still have redundant edges remove the
966 : // duplicates by GrB_mxv and store them as tuples in (SI,SJ,SX)
967 50 : GRB_TRY (GrB_Vector_clear (edge)) ;
968 50 : GRB_TRY (GrB_mxv (edge, mask, tupleMin, minComb, T, I, NULL)) ;
969 50 : GRB_TRY (GrB_Vector_nvals (&num, edge)) ;
970 50 : GRB_TRY (GrB_apply (t, NULL, NULL, get_second, edge, NULL)) ;
971 50 : GRB_TRY (GrB_Vector_extractTuples (NULL, SJ + ntuples, &num, t)) ;
972 50 : GRB_TRY (GrB_apply (w, NULL, NULL, get_first, edge, NULL)) ;
973 :
974 50 : if (weight_type == GrB_INT64)
975 : {
976 28 : GRB_TRY (GrB_Vector_extractTuples_INT64 (
977 : SI + ntuples, ((int64_t *) SX) + ntuples, &num, w)) ;
978 : }
979 : else
980 : {
981 22 : GRB_TRY (GrB_Vector_extractTuples_FP64 (
982 : SI + ntuples, ((double *) SX) + ntuples, &num, w)) ;
983 : }
984 :
985 50 : ntuples += num;
986 :
987 : // path halving until every vertex points on a root
988 : do {
989 152 : GRB_TRY (GxB_Vector_extract_Vector (t, NULL, NULL, f, f, NULL)) ;
990 152 : GRB_TRY (GrB_eWiseMult (mask, NULL, NULL, GrB_NE_UINT64, f, t, NULL)) ;
991 152 : GRB_TRY (GrB_Vector_reduce_BOOL (&diff, NULL, GrB_LOR_MONOID_BOOL, mask, NULL)) ;
992 152 : GrB_Vector temp = f; f = t; t = temp;
993 152 : } while (diff);
994 :
995 : // remove the edges in the same connected component
996 50 : if (weight_type == GrB_INT64)
997 : {
998 28 : GRB_TRY (GrB_Vector_extractTuples (NULL, context_int.parent, &n, f)) ;
999 28 : GRB_TRY (GrB_Matrix_select_UDT (S, NULL, NULL, removeEdge, S, &context_int, NULL)) ;
1000 : }
1001 : else
1002 : {
1003 22 : GRB_TRY (GrB_Vector_extractTuples (NULL, context_fp.parent, &n, f)) ;
1004 22 : GRB_TRY (GrB_Matrix_select_UDT (S, NULL, NULL, removeEdge, S, &context_fp, NULL)) ;
1005 : }
1006 :
1007 50 : GrB_Matrix_nvals (&nvals, S);
1008 : }
1009 :
1010 : // create forest_edges
1011 26 : GRB_TRY (GrB_Matrix_clear (T)) ;
1012 26 : if (weight_type == GrB_INT64)
1013 : {
1014 19 : GRB_TRY (GrB_Matrix_build_INT64 (
1015 : T, SI, SJ, (int64_t *) SX, ntuples, GxB_IGNORE_DUP)) ;
1016 : }
1017 : else
1018 : {
1019 7 : GRB_TRY (GrB_Matrix_build_FP64 (
1020 : T, SI, SJ, (double *) SX, ntuples, GxB_IGNORE_DUP)) ;
1021 : }
1022 :
1023 : //--------------------------------------------------------------------------
1024 : // free workspace and return result
1025 : //--------------------------------------------------------------------------
1026 :
1027 26 : *forest_edges = T;
1028 26 : T = NULL ;
1029 :
1030 26 : if(componentId != NULL)
1031 : {
1032 24 : *componentId = f;
1033 24 : f = NULL;
1034 : }
1035 :
1036 26 : LG_FREE_ALL;
1037 26 : return (GrB_SUCCESS) ;
1038 : #else
1039 : return (GrB_NOT_IMPLEMENTED) ;
1040 : #endif
1041 : }
|
