Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGr_MaxFlow: max flow
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 Darin Peries and Tim Davis, Texas A&M University
15 :
16 : //------------------------------------------------------------------------------
17 :
18 : // LAGr_MaxFlow is a GraphBLAS implementation of the push-relabel algorithm
19 : // of Baumstark et al. [1]
20 : //
21 : // [1] N. Baumstark, G. E. Blelloch, and J. Shun, "Efficient Implementation of
22 : // a Synchronous Parallel Push-Relabel Algorithm." In: Bansal, N., Finocchi, I.
23 : // (eds) Algorithms - ESA 2015. Lecture Notes in Computer Science(), vol 9294.
24 : // Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-662-48350-3 10.
25 :
26 : // [2] D. Peries and T. Davis, "A parallel push-relabel maximum flow algorithm
27 : // in LAGraph and GraphBLAS", IEEE HPEC'25, Sept 2025.
28 :
29 : // TODO: return the (optional) flow matrix can be costly in terms of run time.
30 : // The HPEC'25 results only benchmark the computation of the max flow, f.
31 : // Future work: we plan on revising how the flow matrix is constructed.
32 :
33 : #include <LAGraphX.h>
34 : #include "LG_internal.h"
35 : #include <LAGraph.h>
36 :
37 : //#define DBG
38 : #if LG_SUITESPARSE_GRAPHBLAS_V10
39 :
40 : //------------------------------------------------------------------------------
41 : // LG_augment_maxflow: sum current excess flow into the output flow
42 : //------------------------------------------------------------------------------
43 :
44 : #undef LG_FREE_ALL
45 : #define LG_FREE_ALL ;
46 :
47 362866 : static GrB_Info LG_augment_maxflow
48 : (
49 : double *f, // total maxflow from src to sink
50 : GrB_Vector e, // excess vector, of type double
51 : GrB_Index sink, // sink node
52 : GrB_Vector src_and_sink, // mask vector, with just [src sink]
53 : GrB_Index *n_active, // # of active nodes
54 : char *msg
55 : )
56 : {
57 : // e_sink = e (sink)
58 362866 : double e_sink = 0;
59 362866 : GrB_Info info = GrB_Vector_extractElement(&e_sink, e, sink); //if value at sink
60 362866 : GRB_TRY (info) ;
61 362866 : if (info == GrB_SUCCESS)
62 : {
63 : // e(sink) is present
64 354408 : (*f) += e_sink;
65 : }
66 :
67 : // TODO: what if e is tiny? Do we need a tol parameter,
68 : // and replace all "e > 0" and "r > 0" comparisons with
69 : // (e > tol), throughout the code?
70 :
71 : // e<!struct([src,sink])> = select e where (e > 0)
72 362866 : GRB_TRY(GrB_select(e, src_and_sink, NULL, GrB_VALUEGT_FP64, e, 0, GrB_DESC_RSC));
73 :
74 : // n_active = # of entries in e
75 362866 : GRB_TRY(GrB_Vector_nvals(n_active, e));
76 362866 : return (GrB_SUCCESS) ;
77 : }
78 :
79 : //------------------------------------------------------------------------------
80 : // LG_global_relabel: global relabeling, based on a BFS from the sink node
81 : //------------------------------------------------------------------------------
82 :
83 : #undef LG_FREE_WORK
84 : #define LG_FREE_WORK \
85 : { \
86 : GrB_free(&R_hat); \
87 : GrB_free(&R_hat_transpose); \
88 : LAGraph_Delete(&G2, msg); \
89 : }
90 :
91 : #undef LG_FREE_ALL
92 : #define LG_FREE_ALL LG_FREE_WORK
93 :
94 33904 : static GrB_Info LG_global_relabel
95 : (
96 : // inputs:
97 : GrB_Matrix R, // flow matrix
98 : GrB_Index sink, // sink node
99 : GrB_Vector src_and_sink, // mask vector, with just [src sink]
100 : GrB_UnaryOp GetResidual, // unary op to compute resid=capacity-flow
101 : GrB_BinaryOp global_relabel_accum, // accum for the disconnected assign
102 : GrB_Index relabel_value, // value to relabel the disconnected nodes
103 : // input/output:
104 : GrB_Vector d, // d(i) = height/label of node i
105 : // outputs:
106 : GrB_Vector *lvl, // lvl(i) = distance of node i from sink, if reachable
107 : char *msg
108 : )
109 : {
110 33904 : GrB_Matrix R_hat = NULL, R_hat_transpose = NULL ;
111 33904 : LAGraph_Graph G2 = NULL ;
112 : GrB_Index n ;
113 33904 : GRB_TRY(GrB_Matrix_nrows(&n, R)) ;
114 33904 : GRB_TRY(GrB_Matrix_new(&R_hat_transpose, GrB_FP64, n, n));
115 33904 : GRB_TRY(GrB_Matrix_new(&R_hat, GrB_FP64, n, n));
116 : // R_hat = GetResidual (R), computing the residual of each edge
117 33904 : GRB_TRY(GrB_apply(R_hat, NULL, NULL, GetResidual, R, NULL)) ;
118 : // prune zeros and negative entries from R_hat
119 33904 : GRB_TRY(GrB_select(R_hat, NULL, NULL, GrB_VALUEGT_FP64, R_hat, 0, NULL)) ;
120 : // R_hat_transpose = R_hat'
121 33904 : GRB_TRY(GrB_transpose(R_hat_transpose, NULL, NULL, R_hat, NULL));
122 : // construct G2 and its cached transpose and outdegree
123 33904 : LG_TRY(LAGraph_New(&G2, &R_hat_transpose, LAGraph_ADJACENCY_DIRECTED, msg));
124 33904 : G2->AT = R_hat ;
125 33904 : R_hat = NULL ;
126 33904 : LG_TRY(LAGraph_Cached_OutDegree(G2, msg));
127 : // compute lvl using bfs on G2, starting at sink node
128 33904 : LG_TRY(LAGr_BreadthFirstSearch(lvl, NULL, G2, sink, msg));
129 : // d<!struct([src,sink])> = max(d(i), lvl)
130 33904 : GRB_TRY(GrB_assign(d, src_and_sink, global_relabel_accum, *lvl, GrB_ALL, n, GrB_DESC_SC));
131 : // d<!struct(lvl)> = max(d(i), relabel_value)
132 33904 : GRB_TRY(GrB_assign(d, *lvl, global_relabel_accum, relabel_value, GrB_ALL, n, GrB_DESC_SC));
133 33904 : LG_FREE_WORK ;
134 33904 : return (GrB_SUCCESS) ;
135 : }
136 :
137 : //------------------------------------------------------------------------------
138 :
139 : #define LG_FREE_WORK_EXCEPT_R \
140 : { \
141 : GrB_free(&CompareTuple); \
142 : GrB_free(&e); \
143 : GrB_free(&d); \
144 : GrB_free(&theta); \
145 : GrB_free(&Delta); \
146 : GrB_free(&delta_vec); \
147 : GrB_free(&C); \
148 : GrB_free(&push_vector); \
149 : GrB_free(&pd); \
150 : GrB_free(&src_and_sink); \
151 : GrB_free(&Jvec); \
152 : GrB_free(&Prune); \
153 : GrB_free(&UpdateFlow); \
154 : GrB_free(&Relabel); \
155 : GrB_free(&ResidualFlow); \
156 : GrB_free(&Cxe_IndexMult); \
157 : GrB_free(&Cxe_Mult); \
158 : GrB_free(&Cxe_Semiring); \
159 : GrB_free(&ExtractJ); \
160 : GrB_free(&CreateCompareVec); \
161 : GrB_free(&Rxd_Semiring); \
162 : GrB_free(&Rxd_Add); \
163 : GrB_free(&Rxd_AddMonoid); \
164 : GrB_free(&Rxd_IndexMult); \
165 : GrB_free(&Rxd_Mult); \
166 : GrB_free(&InitForw); \
167 : GrB_free(&InitBack); \
168 : GrB_free(&ResidualForward); \
169 : GrB_free(&ResidualBackward); \
170 : GrB_free(&zero); \
171 : GrB_free(&empty); \
172 : GrB_free(&t); \
173 : GrB_free(&invariant); \
174 : GrB_free(&CheckInvariant); \
175 : GrB_free(&check); \
176 : GrB_free(&ExtractYJ); \
177 : GrB_free(&desc); \
178 : GrB_free(&MakeFlow); \
179 : GrB_free(&GetResidual); \
180 : GrB_free(&lvl) ; \
181 : }
182 :
183 : #undef LG_FREE_WORK
184 : #define LG_FREE_WORK \
185 : { \
186 : LG_FREE_WORK_EXCEPT_R \
187 : GrB_free(&Cxe_Add); \
188 : GrB_free(&Cxe_AddMonoid); \
189 : GrB_free(&ResultTuple); \
190 : GrB_free(&FlowEdge); \
191 : GrB_free(&ExtractMatrixFlow); \
192 : GrB_free(&R); \
193 : }
194 :
195 : #undef LG_FREE_ALL
196 : #define LG_FREE_ALL \
197 : { \
198 : LG_FREE_WORK ; \
199 : GrB_free(flow_mtx) ; \
200 : }
201 :
202 : // helper macro for creating types and operators
203 : #define JIT_STR(f, var) char* var = #f; f
204 :
205 : //casting for unary ops
206 : #define F_UNARY(f) ((void (*)(void *, const void *))f)
207 :
208 : // casting for index binary ops
209 : #define F_INDEX_BINARY(f) ((void (*)(void*, const void*, GrB_Index, GrB_Index, const void *, GrB_Index, GrB_Index, const void *)) f)
210 :
211 : // casting for binary op
212 : #define F_BINARY(f) ((void (*)(void *, const void *, const void *)) f)
213 :
214 : //------------------------------------------------------------------------------
215 : // custom types
216 : //------------------------------------------------------------------------------
217 :
218 : // type of the R matrix: LG_MF_flowEdge (FlowEdge)
219 : JIT_STR(typedef struct{
220 : double capacity; /* original edge weight A(i,j), always positive */
221 : double flow; /* current flow along this edge (i,j); can be negative */
222 : } LG_MF_flowEdge;, FLOWEDGE_STR)
223 :
224 : // type of the push_vector vector for push_vector = R*d: LG_MF_resultTuple64/32 (ResultTuple)
225 : JIT_STR(typedef struct{
226 : double residual; /* residual = capacity - flow for the edge (i,j) */
227 : int64_t d; /* d(j) of the target node j */
228 : int64_t j; /* node id of the target node j */
229 : } LG_MF_resultTuple64;, RESULTTUPLE_STR64)
230 : JIT_STR(typedef struct{
231 : double residual; /* residual = capacity - flow for the edge (i,j) */
232 : int32_t d; /* d(j) of the target node j */
233 : int32_t j; /* node id of the target node j */
234 : } LG_MF_resultTuple32;, RESULTTUPLE_STR32)
235 :
236 : // type of the C matrix and pd vector: LG_MF_compareTuple64/32 (CompareTuple)
237 : JIT_STR(typedef struct{
238 : double residual; /* residual = capacity - flow for the edge (i,j) */
239 : int64_t di; /* d(i) for node i */
240 : int64_t dj; /* d(j) for node j */
241 : int64_t j; /* node id for node j */
242 : } LG_MF_compareTuple64;, COMPARETUPLE_STR64)
243 : JIT_STR(typedef struct{
244 : double residual; /* residual = capacity - flow for the edge (i,j) */
245 : int32_t di; /* d(i) for node i */
246 : int32_t dj; /* d(j) for node j */
247 : int32_t j; /* node id for node j */
248 : int32_t unused; /* to pad the struct to 24 bytes */
249 : } LG_MF_compareTuple32;, COMPARETUPLE_STR32) // 24 bytes: padded
250 :
251 : //------------------------------------------------------------------------------
252 : // unary ops to create R from input adjacency matrix G->A and G->AT
253 : //------------------------------------------------------------------------------
254 :
255 : // unary op for R = ResidualForward (A)
256 19132 : JIT_STR(void LG_MF_ResidualForward(LG_MF_flowEdge *z, const double *y) {
257 : z->capacity = (*y);
258 : z->flow = 0;
259 : }, CRF_STR)
260 :
261 : // unary op for R<!struct(A)> = ResidualBackward (AT)
262 19132 : JIT_STR(void LG_MF_ResidualBackward(LG_MF_flowEdge *z, const double *y) {
263 : z->capacity = 0;
264 : z->flow = 0;
265 : }, CRB_STR)
266 :
267 : //------------------------------------------------------------------------------
268 : // R*d semiring
269 : //------------------------------------------------------------------------------
270 :
271 : // multiplicative operator, z = R(i,j) * d(j), 64-bit case
272 1436 : JIT_STR(void LG_MF_Rxd_Mult64(LG_MF_resultTuple64 *z,
273 : const LG_MF_flowEdge *x, GrB_Index i, GrB_Index j,
274 : const int64_t *y, GrB_Index iy, GrB_Index jy,
275 : const bool* theta) {
276 : double r = x->capacity - x->flow;
277 : if(r > 0){
278 : z->residual = r;
279 : z->d = (*y);
280 : z->j = j;
281 : }
282 : else{
283 : z->residual = 0;
284 : z->d = INT64_MAX;
285 : z->j = -1;
286 : }
287 : }, RXDMULT_STR64)
288 :
289 : // multiplicative operator, z = R(i,j) * d(j), 32-bit case
290 752 : JIT_STR(void LG_MF_Rxd_Mult32(LG_MF_resultTuple32 *z,
291 : const LG_MF_flowEdge *x, GrB_Index i, GrB_Index j,
292 : const int32_t *y, GrB_Index iy, GrB_Index jy,
293 : const bool* theta) {
294 : double r = x->capacity - x->flow;
295 : if(r > 0){
296 : z->residual = r;
297 : z->d = (*y);
298 : z->j = j;
299 : }
300 : else{
301 : z->residual = 0;
302 : z->d = INT32_MAX;
303 : z->j = -1;
304 : }
305 : }, RXDMULT_STR32)
306 :
307 : // additive monoid: z = the best tuple, x or y, 64-bit case
308 1256 : JIT_STR(void LG_MF_Rxd_Add64(LG_MF_resultTuple64 * z,
309 : const LG_MF_resultTuple64 * x,
310 : const LG_MF_resultTuple64 * y) {
311 : if(x->d < y->d){
312 : (*z) = (*x) ;
313 : }
314 : else if(x->d > y->d){
315 : (*z) = (*y) ;
316 : }
317 : else{
318 : if(x->residual > y->residual){
319 : (*z) = (*x) ;
320 : }
321 : else if(x->residual < y->residual){
322 : (*z) = (*y) ;
323 : }
324 : else{
325 : if(x->j > y->j){
326 : (*z) = (*x);
327 : }
328 : else{
329 : (*z) = (*y) ;
330 : }
331 : }
332 : }
333 : }, RXDADD_STR64)
334 :
335 : // additive monoid: z = the best tuple, x or y, 32-bit case
336 464 : JIT_STR(void LG_MF_Rxd_Add32(LG_MF_resultTuple32 * z,
337 : const LG_MF_resultTuple32 * x, const LG_MF_resultTuple32 * y) {
338 : if(x->d < y->d){
339 : (*z) = (*x) ;
340 : }
341 : else if(x->d > y->d){
342 : (*z) = (*y) ;
343 : }
344 : else{
345 : if(x->residual > y->residual){
346 : (*z) = (*x) ;
347 : }
348 : else if(x->residual < y->residual){
349 : (*z) = (*y) ;
350 : }
351 : else{
352 : if(x->j > y->j){
353 : (*z) = (*x);
354 : }
355 : else{
356 : (*z) = (*y) ;
357 : }
358 : }
359 : }
360 : }, RXDADD_STR32)
361 :
362 : //------------------------------------------------------------------------------
363 : // unary ops for delta_vec = ResidualFlow (push_vector)
364 : //------------------------------------------------------------------------------
365 :
366 180 : JIT_STR(void LG_MF_ResidualFlow64(double *z, const LG_MF_resultTuple64 *x)
367 : { (*z) = x->residual; }, RESIDUALFLOW_STR64)
368 :
369 284 : JIT_STR(void LG_MF_ResidualFlow32(double *z, const LG_MF_resultTuple32 *x)
370 : { (*z) = x->residual; }, RESIDUALFLOW_STR32)
371 :
372 : //------------------------------------------------------------------------------
373 : // binary op for R<Delta> = UpdateFlow (R, Delta) using eWiseMult
374 : //------------------------------------------------------------------------------
375 :
376 928 : JIT_STR(void LG_MF_UpdateFlow(LG_MF_flowEdge *z,
377 : const LG_MF_flowEdge *x, const double *y) {
378 : z->capacity = x->capacity;
379 : z->flow = x->flow + (*y);
380 : }, UPDATEFLOW_STR)
381 :
382 : //------------------------------------------------------------------------------
383 : // binary op for d<struct(push_vector)> = Relabel (d, push_vector)
384 : //------------------------------------------------------------------------------
385 :
386 180 : JIT_STR(void LG_MF_Relabel64(int64_t *z,
387 : const int64_t *x, const LG_MF_resultTuple64 *y) {
388 : if((*x) < y->d+1){
389 : (*z) = y->d + 1;
390 : }
391 : else {
392 : (*z) = (*x);
393 : }
394 : }, RELABEL_STR64)
395 :
396 284 : JIT_STR(void LG_MF_Relabel32(int32_t *z,
397 : const int32_t *x, const LG_MF_resultTuple32 *y) {
398 : if((*x) < y->d+1){
399 : (*z) = y->d + 1;
400 : }
401 : else {
402 : (*z) = (*x);
403 : }
404 : }, RELABEL_STR32)
405 :
406 : //------------------------------------------------------------------------------
407 : // unary op for Jvec = ExtractJ (pd), where Jvec(i) = pd(i)->j
408 : //------------------------------------------------------------------------------
409 :
410 180 : JIT_STR(void LG_MF_ExtractJ64(int64_t *z, const LG_MF_compareTuple64 *x) { (*z) = x->j; }, EXTRACTJ_STR64)
411 :
412 288 : JIT_STR(void LG_MF_ExtractJ32(int32_t *z, const LG_MF_compareTuple32 *x) { (*z) = x->j; }, EXTRACTJ_STR32)
413 :
414 : //------------------------------------------------------------------------------
415 : // unary op for Jvec = ExtractYJ(push_vector), where Jvec(i) = push_vector(i)->j
416 : //------------------------------------------------------------------------------
417 :
418 180 : JIT_STR(void LG_MF_ExtractYJ64(int64_t *z, const LG_MF_resultTuple64 *x) { (*z) = x->j; }, EXTRACTYJ_STR64)
419 :
420 284 : JIT_STR(void LG_MF_ExtractYJ32(int32_t *z, const LG_MF_resultTuple32 *x) { (*z) = x->j; }, EXTRACTYJ_STR32)
421 :
422 : //------------------------------------------------------------------------------
423 : // binary op for R(src,:) = InitForw (R (src,:), t')
424 : //------------------------------------------------------------------------------
425 :
426 32 : JIT_STR(void LG_MF_InitForw(LG_MF_flowEdge * z,
427 : const LG_MF_flowEdge * x, const LG_MF_flowEdge * y){
428 : z->capacity = x->capacity;
429 : z->flow = y->flow + x->flow;
430 : }, INITFORW_STR)
431 :
432 : //------------------------------------------------------------------------------
433 : // binary op for R(:,src) = InitBack (R (:,src), t)
434 : //------------------------------------------------------------------------------
435 :
436 32 : JIT_STR(void LG_MF_InitBack(LG_MF_flowEdge * z,
437 : const LG_MF_flowEdge * x, const LG_MF_flowEdge * y){
438 : z->capacity = x->capacity;
439 : z->flow = x->flow - y->flow;
440 : }, INITBACK_STR)
441 :
442 : //------------------------------------------------------------------------------
443 : // push_vector = C*e semiring
444 : //------------------------------------------------------------------------------
445 :
446 : // multiplicative operator, z = C(i,j)*e(j), 64-bit case
447 180 : JIT_STR(void LG_MF_Cxe_Mult64(LG_MF_resultTuple64 * z,
448 : const LG_MF_compareTuple64 * x, GrB_Index i, GrB_Index j,
449 : const double * y, GrB_Index iy, GrB_Index jy,
450 : const bool* theta){
451 : bool j_active = ((*y) > 0) ;
452 : if ((x->di < x->dj-1) /* case a */
453 : || (x->di == x->dj-1 && !j_active) /* case b */
454 : || (x->di == x->dj && (!j_active || (j_active && (i < j)))) /* case c */
455 : || (x->di == x->dj+1)) /* case d */
456 : {
457 : z->residual = x->residual;
458 : z->d = x->dj;
459 : z->j = x->j;
460 : }
461 : else
462 : {
463 : z->residual = 0;
464 : z->d = INT64_MAX;
465 : z->j = -1;
466 : }
467 : }, MXEMULT_STR64)
468 :
469 : // multiplicative operator, z = C(i,j)*e(j), 32-bit case
470 288 : JIT_STR(void LG_MF_Cxe_Mult32(LG_MF_resultTuple32 * z,
471 : const LG_MF_compareTuple32 * x, GrB_Index i, GrB_Index j,
472 : const double * y, GrB_Index iy, GrB_Index jy,
473 : const bool* theta){
474 : bool j_active = ((*y) > 0) ;
475 : if ((x->di < x->dj-1) /* case a */
476 : || (x->di == x->dj-1 && !j_active) /* case b */
477 : || (x->di == x->dj && (!j_active || (j_active && (i < j)))) /* case c */
478 : || (x->di == x->dj+1)) /* case d */
479 : {
480 : z->residual = x->residual;
481 : z->d = x->dj;
482 : z->j = x->j;
483 : }
484 : else
485 : {
486 : z->residual = 0;
487 : z->d = INT32_MAX;
488 : z->j = -1;
489 : }
490 : }, MXEMULT_STR32)
491 :
492 : // Note: the additive monoid is not actually used in the call to GrB_mxv below,
493 : // because any given node only pushes to one neighbor at a time. As a result,
494 : // no reduction is needed in GrB_mxv. The semiring still needs a monoid,
495 : // however.
496 4 : JIT_STR(void LG_MF_Cxe_Add64(LG_MF_resultTuple64 * z,
497 : const LG_MF_resultTuple64 * x, const LG_MF_resultTuple64 * y){
498 : (*z) = (*y) ;
499 : }, MXEADD_STR64)
500 :
501 10 : JIT_STR(void LG_MF_Cxe_Add32(LG_MF_resultTuple32 * z,
502 : const LG_MF_resultTuple32 * x, const LG_MF_resultTuple32 * y){
503 : (*z) = (*y) ;
504 : }, MXEADD_STR32)
505 :
506 : //------------------------------------------------------------------------------
507 : // binary op for pd = CreateCompareVec (push_vector,d) using eWiseMult
508 : //------------------------------------------------------------------------------
509 :
510 180 : JIT_STR(void LG_MF_CreateCompareVec64(LG_MF_compareTuple64 *comp,
511 : const LG_MF_resultTuple64 *res, const int64_t *height) {
512 : comp->di = (*height);
513 : comp->residual = res->residual;
514 : comp->dj = res->d;
515 : comp->j = res->j;
516 : }, CREATECOMPAREVEC_STR64)
517 :
518 288 : JIT_STR(void LG_MF_CreateCompareVec32(LG_MF_compareTuple32 *comp,
519 : const LG_MF_resultTuple32 *res, const int32_t *height) {
520 : comp->di = (*height);
521 : comp->residual = res->residual;
522 : comp->dj = res->d;
523 : comp->j = res->j;
524 : comp->unused = 0 ;
525 : }, CREATECOMPAREVEC_STR32)
526 :
527 : //------------------------------------------------------------------------------
528 : // index unary op to remove empty tuples from push_vector (for which y->j is -1)
529 : //------------------------------------------------------------------------------
530 :
531 686 : JIT_STR(void LG_MF_Prune64(bool * z, const LG_MF_resultTuple64 * x,
532 : GrB_Index ix, GrB_Index jx, const bool * theta){
533 : *z = (x->j != -1) ;
534 : }, PRUNE_STR64)
535 :
536 928 : JIT_STR(void LG_MF_Prune32(bool * z, const LG_MF_resultTuple32 * x,
537 : GrB_Index ix, GrB_Index jx, const bool * theta){
538 : *z = (x->j != -1) ;
539 : }, PRUNE_STR32)
540 :
541 : //------------------------------------------------------------------------------
542 : // unary op for t = MakeFlow (e), where t(i) = (0, e(i))
543 : //------------------------------------------------------------------------------
544 :
545 8080 : JIT_STR(void LG_MF_MakeFlow(LG_MF_flowEdge * flow_edge, const double * flow){
546 : flow_edge->capacity = 0;
547 : flow_edge->flow = (*flow);
548 : }, MAKEFLOW_STR)
549 :
550 : //------------------------------------------------------------------------------
551 : // binary op CheckInvariant to check invariants (debugging only)
552 : //------------------------------------------------------------------------------
553 :
554 : #ifdef DBG
555 : JIT_STR(void LG_MF_CheckInvariant64(bool *z, const int64_t *height,
556 : const LG_MF_resultTuple64 *result) {
557 : (*z) = ((*height) == result->d+1);
558 : }, CHECKINVARIANT_STR64)
559 :
560 : JIT_STR(void LG_MF_CheckInvariant32(bool *z, const int32_t *height,
561 : const LG_MF_resultTuple32 *result) {
562 : (*z) = ((*height) == result->d+1);
563 : }, CHECKINVARIANT_STR32)
564 : #endif
565 :
566 : //------------------------------------------------------------------------------
567 : // binary op for R_hat = GetResidual (R), computing the residual of each edge
568 : //------------------------------------------------------------------------------
569 :
570 14076 : JIT_STR(void LG_MF_GetResidual(double * res, const LG_MF_flowEdge * flow_edge){
571 : (*res) = flow_edge->capacity - flow_edge->flow;
572 : }, GETRESIDUAL_STR)
573 :
574 : //------------------------------------------------------------------------------
575 : // unary op for flow_mtx = ExtractMatrixFlow (R)
576 : //------------------------------------------------------------------------------
577 :
578 5778 : JIT_STR(void LG_MF_ExtractMatrixFlow(double* flow, const LG_MF_flowEdge* edge){
579 : *flow = edge->flow;
580 : }, EMFLOW_STR)
581 :
582 : #endif
583 :
584 : #ifdef DBG
585 : void print_compareVec(const GrB_Vector vec) {
586 : GxB_Iterator iter;
587 : GxB_Iterator_new(&iter);
588 : GrB_Info info = GxB_Vector_Iterator_attach(iter, vec, NULL);
589 : if(info < 0){
590 : printf("error with matrix passed in");
591 : }
592 : info = GxB_Vector_Iterator_seek(iter, 0);
593 : while(info != GxB_EXHAUSTED){
594 : GrB_Index i;
595 : i = GxB_Vector_Iterator_getIndex(iter);
596 : LG_MF_compareTuple32 e;
597 : GxB_Iterator_get_UDT(iter, &e);
598 : printf("(%ld, 0) (di: %d, dj: %d, J: %d, residual: %lf) \n", i, e.di, e.dj, e.j, e.residual);
599 : info = GxB_Vector_Iterator_next(iter);
600 : }
601 : GrB_free(&iter);
602 : }
603 : #endif
604 :
605 :
606 : //------------------------------------------------------------------------------
607 : // LAGraph_MaxFlow
608 : //------------------------------------------------------------------------------
609 :
610 8628 : int LAGr_MaxFlow
611 : (
612 : // output:
613 : double *f, // max flow from src node to sink node
614 : GrB_Matrix *flow_mtx, // optional output flow matrix
615 : // input:
616 : LAGraph_Graph G, // graph to compute maxflow on
617 : GrB_Index src, // source node
618 : GrB_Index sink, // sink node
619 : char *msg
620 : )
621 : {
622 :
623 : #if LG_SUITESPARSE_GRAPHBLAS_V10
624 :
625 : //----------------------------------------------------------------------------
626 : // declare variables
627 : //----------------------------------------------------------------------------
628 :
629 : // types
630 8628 : GrB_Type FlowEdge = NULL ;
631 8628 : GrB_Type ResultTuple = NULL ;
632 8628 : GrB_Type CompareTuple = NULL ;
633 :
634 8628 : GrB_Vector lvl = NULL ;
635 8628 : GrB_UnaryOp GetResidual = NULL ;
636 :
637 : // to create R
638 8628 : GrB_UnaryOp ResidualForward = NULL, ResidualBackward = NULL ;
639 8628 : GrB_Matrix R = NULL ;
640 :
641 : // to initialize R with initial saturated flows
642 8628 : GrB_Vector e = NULL, t = NULL ;
643 8628 : GrB_UnaryOp MakeFlow = NULL ;
644 8628 : GrB_BinaryOp InitForw = NULL, InitBack = NULL ;
645 :
646 : // height/label vector
647 8628 : GrB_Vector d = NULL ;
648 :
649 : // src and sink mask vector and n_active
650 8628 : GrB_Vector src_and_sink = NULL ;
651 8628 : GrB_Index n_active = INT64_MAX ;
652 :
653 : // semiring and vectors for push_vector<struct(e)> = R*d
654 8628 : GrB_Vector push_vector = NULL ;
655 8628 : GrB_IndexUnaryOp Prune = NULL ;
656 8628 : GxB_IndexBinaryOp Rxd_IndexMult = NULL ;
657 8628 : GrB_BinaryOp Rxd_Add = NULL, Rxd_Mult = NULL ;
658 8628 : GrB_Monoid Rxd_AddMonoid = NULL ;
659 8628 : GrB_Semiring Rxd_Semiring = NULL ;
660 8628 : GrB_Scalar theta = NULL ;
661 :
662 : // binary op and pd
663 8628 : GrB_Vector pd = NULL ;
664 8628 : GrB_BinaryOp CreateCompareVec = NULL ;
665 :
666 : // utility vectors, Matrix, and ops for mapping
667 8628 : GrB_Matrix C = NULL ; // matrix of candidate pushes
668 8628 : GrB_Vector Jvec = NULL ;
669 8628 : GrB_UnaryOp ExtractJ = NULL, ExtractYJ = NULL ;
670 :
671 : // C*e semiring
672 8628 : GrB_Semiring Cxe_Semiring = NULL ;
673 8628 : GrB_Monoid Cxe_AddMonoid = NULL ;
674 8628 : GrB_BinaryOp Cxe_Add = NULL, Cxe_Mult = NULL ;
675 8628 : GxB_IndexBinaryOp Cxe_IndexMult = NULL ;
676 :
677 : // to extract the residual flow
678 8628 : GrB_UnaryOp ResidualFlow = NULL ;
679 :
680 : // to extract the final flows for the flow matrix
681 8628 : GrB_UnaryOp ExtractMatrixFlow = NULL ;
682 :
683 : // Delta structures
684 8628 : GrB_Vector delta_vec = NULL ;
685 8628 : GrB_Matrix Delta = NULL ;
686 :
687 : // update height
688 8628 : GrB_BinaryOp Relabel = NULL ;
689 :
690 : // update R structure
691 8628 : GrB_BinaryOp UpdateFlow = NULL ;
692 :
693 : // scalars
694 8628 : GrB_Scalar zero = NULL ;
695 8628 : GrB_Scalar empty = NULL ;
696 :
697 : // invariant (for debugging only)
698 8628 : GrB_Vector invariant = NULL ;
699 8628 : GrB_BinaryOp CheckInvariant = NULL ;
700 8628 : GrB_Scalar check = NULL ;
701 : bool check_raw;
702 :
703 : // descriptor for matrix building
704 8628 : GrB_Descriptor desc = NULL ;
705 :
706 : //----------------------------------------------------------------------------
707 : // check inputs
708 : //----------------------------------------------------------------------------
709 :
710 8628 : if (flow_mtx != NULL)
711 : {
712 14 : (*flow_mtx) = NULL ;
713 : }
714 8628 : LG_TRY(LAGraph_CheckGraph(G, msg));
715 8628 : LG_ASSERT (f != NULL, GrB_NULL_POINTER) ;
716 8628 : (*f) = 0;
717 : GrB_Index nrows, n;
718 8628 : GRB_TRY(GrB_Matrix_ncols(&n, G->A));
719 8628 : GRB_TRY(GrB_Matrix_nrows(&nrows, G->A));
720 8628 : LG_ASSERT_MSG(nrows == n, GrB_INVALID_VALUE, "Matrix must be square");
721 8628 : LG_ASSERT_MSG(src < n && src >= 0 && sink < n && sink >= 0,
722 : GrB_INVALID_VALUE, "src and sink must be a value between [0, n)");
723 8628 : LG_ASSERT_MSG(G->emin > 0, GrB_INVALID_VALUE,
724 : "the edge weights (capacities) must be greater than 0");
725 :
726 : //get adjacency matrix and its transpose
727 8628 : GrB_Matrix A = G->A;
728 8628 : GrB_Matrix AT = NULL ;
729 8628 : if (G->kind == LAGraph_ADJACENCY_UNDIRECTED)
730 : {
731 : // G is undirected, so A and AT are the same
732 4 : AT = G->A ;
733 : }
734 : else
735 : {
736 : // G is directed; get G->AT, which must be present
737 8624 : AT = G->AT ;
738 8624 : LG_ASSERT_MSG (AT != NULL, LAGRAPH_NOT_CACHED, "G->AT is required") ;
739 : }
740 :
741 : //----------------------------------------------------------------------------
742 : // create types, operators, matrices, and vectors
743 : //----------------------------------------------------------------------------
744 :
745 : // create types for computation
746 8628 : GRB_TRY(GxB_Type_new(&FlowEdge, sizeof(LG_MF_flowEdge), "LG_MF_flowEdge", FLOWEDGE_STR));
747 :
748 8628 : GRB_TRY(GxB_UnaryOp_new(&GetResidual, F_UNARY(LG_MF_GetResidual), GrB_FP64, FlowEdge,
749 : "LG_MF_GetResidual", GETRESIDUAL_STR));
750 :
751 : #ifdef DBG
752 : GRB_TRY(GrB_Scalar_new(&check, GrB_BOOL));
753 : GRB_TRY(GrB_Scalar_setElement(check, false));
754 : GRB_TRY(GrB_Vector_new(&invariant, GrB_BOOL, n));
755 : #endif
756 :
757 : // ops create R from A
758 8628 : GRB_TRY(GxB_UnaryOp_new(&ResidualForward,
759 : F_UNARY(LG_MF_ResidualForward), FlowEdge , GrB_FP64,
760 : "LG_MF_ResidualForward", CRF_STR));
761 8628 : GRB_TRY(GxB_UnaryOp_new(&ResidualBackward,
762 : F_UNARY(LG_MF_ResidualBackward), FlowEdge , GrB_FP64,
763 : "LG_MF_ResidualBackward", CRB_STR));
764 :
765 : // ops to initialize R with initial saturated flows from the source node
766 8628 : GRB_TRY(GxB_BinaryOp_new(&InitForw,
767 : F_BINARY(LG_MF_InitForw), FlowEdge, FlowEdge, FlowEdge,
768 : "LG_MF_InitForw", INITFORW_STR));
769 8628 : GRB_TRY(GxB_BinaryOp_new(&InitBack,
770 : F_BINARY(LG_MF_InitBack), FlowEdge, FlowEdge, FlowEdge,
771 : "LG_MF_InitBack", INITBACK_STR));
772 8628 : GRB_TRY(GxB_UnaryOp_new(&MakeFlow, F_UNARY(LG_MF_MakeFlow), FlowEdge, GrB_FP64,
773 : "LG_MF_MakeFlow", MAKEFLOW_STR));
774 :
775 : // construct [src,sink] mask
776 8628 : GRB_TRY(GrB_Vector_new(&src_and_sink, GrB_BOOL, n));
777 8628 : GRB_TRY (GrB_Vector_setElement (src_and_sink, true, sink)) ;
778 8628 : GRB_TRY (GrB_Vector_setElement (src_and_sink, true, src)) ;
779 :
780 : // create delta vector and Delta matrix
781 8628 : GRB_TRY(GrB_Matrix_new(&Delta, GrB_FP64, n, n));
782 8628 : GRB_TRY(GrB_Vector_new(&delta_vec, GrB_FP64, n));
783 :
784 : // operator to update R structure
785 8628 : GRB_TRY(GxB_BinaryOp_new(&UpdateFlow,
786 : F_BINARY(LG_MF_UpdateFlow), FlowEdge, FlowEdge, GrB_FP64,
787 : "LG_MF_UpdateFlow", UPDATEFLOW_STR));
788 :
789 : // create scalars
790 8628 : GRB_TRY(GrB_Scalar_new(&zero, GrB_FP64));
791 8628 : GRB_TRY(GrB_Scalar_setElement(zero, 0));
792 8628 : GRB_TRY(GrB_Scalar_new (&empty, GrB_FP64)) ;
793 8628 : GRB_TRY(GrB_Scalar_new(&theta, GrB_BOOL)); // unused placeholder
794 8628 : GRB_TRY(GrB_Scalar_setElement_BOOL(theta, false));
795 :
796 : //----------------------------------------------------------------------------
797 : // determine the integer type to use for the problem
798 : //----------------------------------------------------------------------------
799 :
800 8628 : GrB_Type Integer_Type = NULL ;
801 :
802 : //accum operator for the global relabel
803 8628 : GrB_BinaryOp global_relabel_accum = NULL ;
804 :
805 : #ifdef COVERAGE
806 : // Just for test coverage, use 64-bit ints for n > 100. Do not use this
807 : // rule in production!
808 : #define NBIG 100
809 : #else
810 : // For production use: 64-bit integers if n > 2^31
811 : #define NBIG INT32_MAX
812 : #endif
813 8628 : if (n > NBIG){
814 :
815 : //--------------------------------------------------------------------------
816 : // use 64-bit integers
817 : //--------------------------------------------------------------------------
818 :
819 8 : Integer_Type = GrB_INT64 ;
820 :
821 : // use the 64 bit max operator
822 8 : global_relabel_accum = GrB_MAX_INT64 ;
823 :
824 : // create types for computation
825 8 : GRB_TRY(GxB_Type_new(&ResultTuple, sizeof(LG_MF_resultTuple64),
826 : "LG_MF_resultTuple64", RESULTTUPLE_STR64));
827 8 : GRB_TRY(GxB_Type_new(&CompareTuple, sizeof(LG_MF_compareTuple64),
828 : "LG_MF_compareTuple64", COMPARETUPLE_STR64));
829 :
830 : // invariant check
831 : #ifdef DBG
832 : GRB_TRY(GxB_BinaryOp_new(&CheckInvariant,
833 : F_BINARY(LG_MF_CheckInvariant64), GrB_BOOL, GrB_INT64, ResultTuple,
834 : "LG_MF_CheckInvariant64", CHECKINVARIANT_STR64));
835 : #endif
836 :
837 8 : GRB_TRY(GxB_UnaryOp_new(&ResidualFlow,
838 : F_UNARY(LG_MF_ResidualFlow64), GrB_FP64, ResultTuple,
839 : "LG_MF_ResidualFlow64", RESIDUALFLOW_STR64));
840 :
841 : // create ops for R*d semiring
842 :
843 8 : GRB_TRY(GxB_IndexBinaryOp_new(&Rxd_IndexMult,
844 : F_INDEX_BINARY(LG_MF_Rxd_Mult64), ResultTuple, FlowEdge, GrB_INT64, GrB_BOOL,
845 : "LG_MF_Rxd_Mult64", RXDMULT_STR64));
846 8 : GRB_TRY(GxB_BinaryOp_new_IndexOp(&Rxd_Mult, Rxd_IndexMult, theta));
847 8 : GRB_TRY(GxB_BinaryOp_new(&Rxd_Add,
848 : F_BINARY(LG_MF_Rxd_Add64), ResultTuple, ResultTuple, ResultTuple,
849 : "LG_MF_Rxd_Add64", RXDADD_STR64));
850 8 : LG_MF_resultTuple64 id = {.d = INT64_MAX, .j = -1, .residual = 0};
851 8 : GRB_TRY(GrB_Monoid_new_UDT(&Rxd_AddMonoid, Rxd_Add, &id));
852 :
853 : // create binary op for pd
854 8 : GRB_TRY(GxB_BinaryOp_new(&CreateCompareVec,
855 : F_BINARY(LG_MF_CreateCompareVec64), CompareTuple, ResultTuple, GrB_INT64,
856 : "LG_MF_CreateCompareVec64", CREATECOMPAREVEC_STR64));
857 :
858 : // create op to prune empty tuples
859 8 : GRB_TRY(GxB_IndexUnaryOp_new(&Prune,
860 : (GxB_index_unary_function) LG_MF_Prune64, GrB_BOOL, ResultTuple, GrB_BOOL,
861 : "LG_MF_Prune64", PRUNE_STR64));
862 :
863 : // create ops for mapping
864 8 : GRB_TRY(GxB_UnaryOp_new(&ExtractJ,
865 : F_UNARY(LG_MF_ExtractJ64), GrB_INT64, CompareTuple,
866 : "LG_MF_ExtractJ64", EXTRACTJ_STR64));
867 8 : GRB_TRY(GxB_UnaryOp_new(&ExtractYJ,
868 : F_UNARY(LG_MF_ExtractYJ64), GrB_INT64, ResultTuple,
869 : "LG_MF_ExtractYJ64", EXTRACTYJ_STR64));
870 :
871 : // create ops for C*e semiring
872 8 : GRB_TRY(GxB_IndexBinaryOp_new(&Cxe_IndexMult,
873 : F_INDEX_BINARY(LG_MF_Cxe_Mult64), ResultTuple, CompareTuple, GrB_FP64, GrB_BOOL,
874 : "LG_MF_Cxe_Mult64", MXEMULT_STR64));
875 8 : GRB_TRY(GxB_BinaryOp_new_IndexOp(&Cxe_Mult, Cxe_IndexMult, theta));
876 8 : GRB_TRY(GxB_BinaryOp_new(&Cxe_Add,
877 : F_BINARY(LG_MF_Cxe_Add64), ResultTuple, ResultTuple, ResultTuple,
878 : "LG_MF_Cxe_Add64", MXEADD_STR64));
879 8 : GRB_TRY(GrB_Monoid_new_UDT(&Cxe_AddMonoid, Cxe_Add, &id));
880 :
881 : // update height binary op
882 8 : GRB_TRY(GxB_BinaryOp_new(&Relabel,
883 : F_BINARY(LG_MF_Relabel64), GrB_INT64, GrB_INT64, ResultTuple,
884 : "LG_MF_Relabel64", RELABEL_STR64));
885 :
886 : }else{
887 :
888 : //--------------------------------------------------------------------------
889 : // use 32-bit integers
890 : //--------------------------------------------------------------------------
891 :
892 8620 : Integer_Type = GrB_INT32 ;
893 :
894 : // use 32 bit max op
895 8620 : global_relabel_accum = GrB_MAX_INT32 ;
896 :
897 : // create types for computation
898 8620 : GRB_TRY(GxB_Type_new(&ResultTuple, sizeof(LG_MF_resultTuple32),
899 : "LG_MF_resultTuple32", RESULTTUPLE_STR32));
900 8620 : GRB_TRY(GxB_Type_new(&CompareTuple, sizeof(LG_MF_compareTuple32),
901 : "LG_MF_compareTuple32", COMPARETUPLE_STR32));
902 :
903 : // invariant check
904 : #ifdef DBG
905 : GRB_TRY(GxB_BinaryOp_new(&CheckInvariant,
906 : F_BINARY(LG_MF_CheckInvariant32), GrB_BOOL, GrB_INT32, ResultTuple,
907 : "LG_MF_CheckInvariant32", CHECKINVARIANT_STR32));
908 : #endif
909 :
910 8620 : GRB_TRY(GxB_UnaryOp_new(&ResidualFlow,
911 : F_UNARY(LG_MF_ResidualFlow32), GrB_FP64, ResultTuple,
912 : "LG_MF_ResidualFlow32", RESIDUALFLOW_STR32));
913 :
914 : // create ops for R*d semiring
915 8620 : GRB_TRY(GxB_IndexBinaryOp_new(&Rxd_IndexMult,
916 : F_INDEX_BINARY(LG_MF_Rxd_Mult32), ResultTuple, FlowEdge, GrB_INT32, GrB_BOOL,
917 : "LG_MF_Rxd_Mult32", RXDMULT_STR32));
918 8620 : GRB_TRY(GxB_BinaryOp_new_IndexOp(&Rxd_Mult, Rxd_IndexMult, theta));
919 8620 : GRB_TRY(GxB_BinaryOp_new(&Rxd_Add,
920 : F_BINARY(LG_MF_Rxd_Add32), ResultTuple, ResultTuple, ResultTuple,
921 : "LG_MF_Rxd_Add32", RXDADD_STR32));
922 8620 : LG_MF_resultTuple32 id = {.d = INT32_MAX, .j = -1, .residual = 0};
923 8620 : GRB_TRY(GrB_Monoid_new_UDT(&Rxd_AddMonoid, Rxd_Add, &id));
924 :
925 : // create binary op for pd
926 8620 : GRB_TRY(GxB_BinaryOp_new(&CreateCompareVec,
927 : F_BINARY(LG_MF_CreateCompareVec32), CompareTuple, ResultTuple, GrB_INT32,
928 : "LG_MF_CreateCompareVec32", CREATECOMPAREVEC_STR32));
929 :
930 : // create op to prune empty tuples
931 8620 : GRB_TRY(GxB_IndexUnaryOp_new(&Prune,
932 : (GxB_index_unary_function) LG_MF_Prune32, GrB_BOOL, ResultTuple, GrB_BOOL,
933 : "LG_MF_Prune32", PRUNE_STR32));
934 :
935 : // create ops for mapping
936 8620 : GRB_TRY(GxB_UnaryOp_new(&ExtractJ,
937 : F_UNARY(LG_MF_ExtractJ32), GrB_INT32, CompareTuple,
938 : "LG_MF_ExtractJ32", EXTRACTJ_STR32));
939 8620 : GRB_TRY(GxB_UnaryOp_new(&ExtractYJ,
940 : F_UNARY(LG_MF_ExtractYJ32), GrB_INT32, ResultTuple,
941 : "LG_MF_ExtractYJ32", EXTRACTYJ_STR32));
942 :
943 : // create ops for C*e semiring
944 8620 : GRB_TRY(GxB_IndexBinaryOp_new(&Cxe_IndexMult,
945 : F_INDEX_BINARY(LG_MF_Cxe_Mult32), ResultTuple, CompareTuple, GrB_FP64, GrB_BOOL,
946 : "LG_MF_Cxe_Mult32", MXEMULT_STR32));
947 8620 : GRB_TRY(GxB_BinaryOp_new_IndexOp(&Cxe_Mult, Cxe_IndexMult, theta));
948 8620 : GRB_TRY(GxB_BinaryOp_new(&Cxe_Add,
949 : F_BINARY(LG_MF_Cxe_Add32), ResultTuple, ResultTuple, ResultTuple,
950 : "LG_MF_Cxe_Add32", MXEADD_STR32));
951 8620 : GRB_TRY(GrB_Monoid_new_UDT(&Cxe_AddMonoid, Cxe_Add, &id));
952 :
953 : // update height binary op
954 8620 : GRB_TRY(GxB_BinaryOp_new(&Relabel,
955 : F_BINARY(LG_MF_Relabel32), GrB_INT32, GrB_INT32, ResultTuple,
956 : "LG_MF_Relabel32", RELABEL_STR32));
957 : }
958 :
959 : //----------------------------------------------------------------------------
960 : // create remaining vectors, matrices, descriptor, and semirings
961 : //----------------------------------------------------------------------------
962 :
963 8628 : GRB_TRY(GrB_Matrix_new(&C, CompareTuple, n,n));
964 8628 : GRB_TRY(GrB_Vector_new(&Jvec, Integer_Type, n));
965 8628 : GRB_TRY(GrB_Vector_new(&pd, CompareTuple, n));
966 8628 : GRB_TRY(GrB_Vector_new(&push_vector, ResultTuple, n));
967 :
968 8628 : GRB_TRY(GrB_Semiring_new(&Rxd_Semiring, Rxd_AddMonoid, Rxd_Mult));
969 8628 : GRB_TRY(GrB_Semiring_new(&Cxe_Semiring, Cxe_AddMonoid, Cxe_Mult));
970 :
971 : // create descriptor for building the C and Delta matrices
972 8628 : GRB_TRY(GrB_Descriptor_new(&desc));
973 8628 : GRB_TRY(GrB_set(desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST));
974 :
975 : // create and init d vector
976 8628 : GRB_TRY(GrB_Vector_new(&d, Integer_Type, n));
977 8628 : GRB_TRY(GrB_assign(d, NULL, NULL, 0, GrB_ALL, n, NULL));
978 8628 : GRB_TRY(GrB_assign(d, NULL, NULL, n, &src, 1, NULL));
979 :
980 : // create R, with no flow
981 8628 : GRB_TRY(GrB_Matrix_new(&R, FlowEdge, n, n));
982 : // R = ResidualForward (A)
983 8628 : GRB_TRY(GrB_apply(R, NULL, NULL, ResidualForward, A, NULL));
984 : // R<!struct(A)> = ResidualBackward (AT)
985 8628 : GRB_TRY(GrB_apply(R, A, NULL, ResidualBackward, AT, GrB_DESC_SC));
986 :
987 : // initial global relabeling
988 : // relabel to 2*n to prevent any flow from going to the
989 : // disconnected nodes.
990 8628 : GrB_Index relabel_value = 2*n ;
991 8628 : LG_TRY (LG_global_relabel (R, sink, src_and_sink, GetResidual, global_relabel_accum, relabel_value, d, &lvl, msg)) ;
992 :
993 : // reset to n
994 8628 : relabel_value = n ;
995 :
996 : // create excess vector e and initial flows from the src to its neighbors
997 : // e<struct(lvl)> = A (src,:)
998 8628 : GRB_TRY(GrB_Vector_new(&e, GrB_FP64, n));
999 8628 : GRB_TRY(GrB_extract(e, lvl, NULL, A, GrB_ALL, n, src, GrB_DESC_ST0));
1000 8628 : GrB_free(&lvl);
1001 : // t = MakeFlow (e), where t(i) = (0, e(i))
1002 8628 : GRB_TRY(GrB_Vector_new(&t, FlowEdge, n));
1003 8628 : GRB_TRY(GrB_apply(t, NULL, NULL, MakeFlow, e, NULL));
1004 : // R(src,:) = InitForw (R (src,:), t')
1005 8628 : GRB_TRY(GrB_assign(R, NULL, InitForw, t, src, GrB_ALL, n, NULL));
1006 : // R(:,src) = InitBack (R (:,src), t)
1007 8628 : GRB_TRY(GrB_assign(R, NULL, InitBack, t, GrB_ALL, n, src, NULL));
1008 8628 : GrB_free(&t) ;
1009 :
1010 : // augment the maxflow with the initial flows from the src to its neighbors
1011 8628 : LG_TRY (LG_augment_maxflow (f, e, sink, src_and_sink, &n_active, msg)) ;
1012 :
1013 : //----------------------------------------------------------------------------
1014 : // compute the max flow
1015 : //----------------------------------------------------------------------------
1016 :
1017 362866 : for (int64_t iter = 0 ; n_active > 0 ; iter++)
1018 : {
1019 : #ifdef DBG
1020 : printf ("iter: %ld, n_active %ld\n", iter, n_active) ;
1021 : #endif
1022 : //--------------------------------------------------------------------------
1023 : // Part 1: global relabeling
1024 : //--------------------------------------------------------------------------
1025 :
1026 354282 : if ((iter > 0) && (iter % 12 == 0))
1027 : {
1028 : #ifdef DBG
1029 : printf ("relabel at : %ld\n", iter) ;
1030 : #endif
1031 25276 : LG_TRY (LG_global_relabel (R, sink, src_and_sink, GetResidual, global_relabel_accum, relabel_value, d, &lvl, msg)) ;
1032 25276 : if(flow_mtx == NULL){
1033 : // delete nodes in e that cannot be reached from the sink
1034 : // e<!struct(lvl)> = empty scalar
1035 25262 : GrB_assign (e, lvl, NULL, empty, GrB_ALL, n, GrB_DESC_SC) ;
1036 25262 : GRB_TRY(GrB_Vector_nvals(&n_active, e));
1037 25262 : if(n_active == 0) break;
1038 : }
1039 25232 : GrB_free(&lvl);
1040 : }
1041 :
1042 : //--------------------------------------------------------------------------
1043 : // Part 2: deciding where to push
1044 : //--------------------------------------------------------------------------
1045 :
1046 : // push_vector<struct(e),replace> = R*d using the Rxd_Semiring
1047 354238 : GRB_TRY(GrB_mxv(push_vector, e, NULL, Rxd_Semiring, R, d, GrB_DESC_RS));
1048 :
1049 : // remove empty tuples (0,inf,-1) from push_vector
1050 354238 : GRB_TRY(GrB_select(push_vector, NULL, NULL, Prune, push_vector, 0, NULL));
1051 :
1052 : //--------------------------------------------------------------------------
1053 : // Part 3: verifying the pushes
1054 : //--------------------------------------------------------------------------
1055 :
1056 : // create C matrix (Candidate pushes) from pattern and values of pd
1057 : // pd = CreateCompareVec (push_vector,d) using eWiseMult
1058 354238 : GRB_TRY(GrB_eWiseMult(pd, NULL, NULL, CreateCompareVec, push_vector, d, NULL));
1059 :
1060 : #ifdef DBG
1061 : print_compareVec(pd);
1062 : GxB_print(d, 3);
1063 : GxB_print(e, 3);
1064 : #endif
1065 :
1066 : // Jvec = ExtractJ (pd), where Jvec(i) = pd(i)->j
1067 354238 : GRB_TRY(GrB_apply(Jvec, NULL, NULL, ExtractJ, pd, NULL));
1068 354238 : GRB_TRY(GrB_Matrix_clear(C));
1069 354238 : GRB_TRY(GrB_Matrix_build(C, pd, Jvec, pd, GxB_IGNORE_DUP, desc));
1070 354238 : GRB_TRY(GrB_Vector_clear(pd));
1071 354238 : GRB_TRY(GrB_Vector_clear(Jvec));
1072 :
1073 : // make e dense for C computation
1074 : // TODO: consider keeping e in bitmap/full format only,
1075 : // or always full with e(i)=0 denoting a non-active node.
1076 354238 : GRB_TRY(GrB_assign(e, e, NULL, 0, GrB_ALL, n, GrB_DESC_SC));
1077 :
1078 : // push_vector = C*e using the Cxe_Semiring
1079 354238 : GRB_TRY(GrB_mxv(push_vector, NULL, NULL, Cxe_Semiring, C, e, NULL));
1080 354238 : GRB_TRY(GrB_Matrix_clear(C));
1081 :
1082 : // remove empty tuples (0,inf,-1) from push_vector
1083 354238 : GRB_TRY(GrB_select(push_vector, NULL, NULL, Prune, push_vector, -1, NULL));
1084 :
1085 : // relabel, updating the height/label vector d
1086 : // d<struct(push_vector)> = Relabel (d, push_vector) using eWiseMult
1087 354238 : GRB_TRY(GrB_eWiseMult(d, push_vector, NULL, Relabel, d, push_vector, GrB_DESC_S));
1088 :
1089 : #ifdef DBG
1090 : // assert invariant for all labels
1091 : GRB_TRY(GrB_eWiseMult(invariant, push_vector, NULL, CheckInvariant, d, push_vector, GrB_DESC_RS));
1092 : GRB_TRY(GrB_reduce(check, NULL, GrB_LAND_MONOID_BOOL, invariant, NULL));
1093 : GRB_TRY(GrB_Scalar_extractElement(&check_raw, check));
1094 : ASSERT(check_raw == true);
1095 : #endif
1096 :
1097 : //--------------------------------------------------------------------------
1098 : // Part 4: executing the pushes
1099 : //--------------------------------------------------------------------------
1100 :
1101 : // extract residual flows from push_vector
1102 : // delta_vec = ResidualFlow (push_vector), obtaining just the residual flows
1103 354238 : GRB_TRY(GrB_apply(delta_vec, NULL, NULL, ResidualFlow, push_vector, NULL));
1104 :
1105 : // delta_vec = min (delta_vec, e), where e is dense
1106 354238 : GRB_TRY(GrB_eWiseMult(delta_vec, NULL, NULL, GrB_MIN_FP64, delta_vec, e, NULL));
1107 :
1108 : // create the Delta matrix from delta_vec and push_vector
1109 : // note that delta_vec has the same structure as push_vector
1110 : // Jvec = ExtractYJ (push_vector), where Jvec(i) = push_vector(i)->j
1111 : // if Jvec has no values, then there is no possible
1112 : // candidates to push to, so the algorithm terminates
1113 354238 : GRB_TRY(GrB_apply(Jvec, NULL, NULL, ExtractYJ, push_vector, NULL));
1114 : GrB_Index J_n;
1115 354238 : GRB_TRY(GrB_Vector_nvals(&J_n, Jvec));
1116 354238 : if(J_n == 0) break;
1117 354238 : GRB_TRY(GrB_Matrix_clear(Delta));
1118 354238 : GRB_TRY(GrB_Matrix_build(Delta, delta_vec, Jvec, delta_vec, GxB_IGNORE_DUP, desc));
1119 354238 : GRB_TRY(GrB_Vector_clear(Jvec));
1120 :
1121 : // make Delta anti-symmetric
1122 : // Delta = (Delta - Delta')
1123 354238 : GRB_TRY(GxB_eWiseUnion(Delta, NULL, NULL, GrB_MINUS_FP64, Delta, zero, Delta, zero, GrB_DESC_T1));
1124 :
1125 : // update R
1126 : // R<Delta> = UpdateFlow (R, Delta) using eWiseMult
1127 354238 : GRB_TRY(GrB_eWiseMult(R, Delta, NULL, UpdateFlow, R, Delta, GrB_DESC_S));
1128 :
1129 : // reduce Delta to delta_vec
1130 : // delta_vec = sum (Delta), summing up each row of Delta
1131 354238 : GRB_TRY(GrB_reduce(delta_vec, NULL, NULL, GrB_PLUS_MONOID_FP64, Delta, GrB_DESC_T0));
1132 354238 : GRB_TRY(GrB_Matrix_clear(Delta));
1133 :
1134 : // add delta_vec to e
1135 : // e<struct(delta_vec)> += delta_vec
1136 354238 : GRB_TRY(GrB_assign(e, delta_vec, GrB_PLUS_FP64, delta_vec, GrB_ALL, n, GrB_DESC_S));
1137 :
1138 : // augment maxflow for all active nodes
1139 354238 : LG_TRY (LG_augment_maxflow (f, e, sink, src_and_sink, &n_active, msg)) ;
1140 :
1141 : }
1142 :
1143 :
1144 : //----------------------------------------------------------------------------
1145 : // optionally construct the output flow matrix, if requested
1146 : //----------------------------------------------------------------------------
1147 :
1148 8628 : if (flow_mtx != NULL)
1149 : {
1150 : // free all workspace except R
1151 14 : LG_FREE_WORK_EXCEPT_R ;
1152 : // create the ExtractMatrixFlow op to compute the flow matrix
1153 14 : GRB_TRY(GxB_UnaryOp_new(&ExtractMatrixFlow,
1154 : F_UNARY(LG_MF_ExtractMatrixFlow), GrB_FP64, FlowEdge,
1155 : "LG_MF_ExtractMatrixFlow", EMFLOW_STR));
1156 : // flow_mtx = ExtractMatrixFlow (R)
1157 14 : GRB_TRY(GrB_Matrix_new(flow_mtx, GrB_FP64, n, n));
1158 14 : GRB_TRY(GrB_apply(*flow_mtx, NULL, NULL, ExtractMatrixFlow, R, NULL));
1159 : // delete any zero or negative flows from the flow_mtx
1160 14 : GRB_TRY(GrB_select(*flow_mtx, NULL, NULL, GrB_VALUEGT_FP64, *flow_mtx, 0, NULL));
1161 : }
1162 :
1163 : //----------------------------------------------------------------------------
1164 : // for test coverage only
1165 : //----------------------------------------------------------------------------
1166 :
1167 : #ifdef COVERAGE
1168 : // The Cxe_Add operator is not tested via the call to GrB_mxv with the
1169 : // Cxe_Semiring above, so test it via the Cxe_AddMonoid.
1170 8628 : GrB_free(&push_vector);
1171 8628 : GRB_TRY(GrB_Vector_new(&push_vector, ResultTuple, 3));
1172 8628 : if (n > NBIG)
1173 : {
1174 8 : LG_MF_resultTuple64 a = {.d = 1, .j = 2, .residual = 3};
1175 8 : LG_MF_resultTuple64 b = {.d = 4, .j = 5, .residual = 6};
1176 8 : GRB_TRY (GrB_Vector_setElement_UDT (push_vector, (void *) &a, 0)) ;
1177 8 : GRB_TRY (GrB_Vector_setElement_UDT (push_vector, (void *) &b, 0)) ;
1178 8 : LG_MF_resultTuple64 c = {.d = 0, .j = 0, .residual = 0};
1179 8 : GRB_TRY (GrB_Vector_reduce_UDT ((void *) &c, NULL, Cxe_AddMonoid, push_vector, NULL)) ;
1180 8 : LG_ASSERT ((c.residual == 6 && c.j == 5 && c.d == 4), GrB_PANIC) ;
1181 : }
1182 : else
1183 : {
1184 8620 : LG_MF_resultTuple32 a = {.d = 1, .j = 2, .residual = 3};
1185 8620 : LG_MF_resultTuple32 b = {.d = 4, .j = 5, .residual = 6};
1186 8620 : GRB_TRY (GrB_Vector_setElement_UDT (push_vector, (void *) &a, 0)) ;
1187 8620 : GRB_TRY (GrB_Vector_setElement_UDT (push_vector, (void *) &b, 0)) ;
1188 8620 : LG_MF_resultTuple32 c = {.d = 0, .j = 0, .residual = 0};
1189 8620 : GRB_TRY (GrB_Vector_reduce_UDT ((void *) &c, NULL, Cxe_AddMonoid, push_vector, NULL)) ;
1190 8620 : LG_ASSERT ((c.residual == 6 && c.j == 5 && c.d == 4), GrB_PANIC) ;
1191 : }
1192 : #endif
1193 :
1194 : //----------------------------------------------------------------------------
1195 : // free workspace and return result
1196 : //----------------------------------------------------------------------------
1197 :
1198 8628 : LG_FREE_WORK ;
1199 8628 : return GrB_SUCCESS;
1200 : #else
1201 : return GrB_NOT_IMPLEMENTED ;
1202 : #endif
1203 : }
|
