Skip to content

Commit 5d17fa8

Browse files
DrTimothyAldenDavisdperry17
DrTimothyAldenDavis
authored and
dperry17
committed
minor cleanup. Add FLOP tags and a TODO for e > tol
1 parent 43c0fcf commit 5d17fa8

File tree

1 file changed

+84
-53
lines changed

1 file changed

+84
-53
lines changed

‎experimental/algorithm/LAGr_MaxFlow.c‎

Lines changed: 84 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -21,29 +21,48 @@
2121

2222
static GrB_Info LG_augment_maxflow
2323
(
24-
double *f, // maxflow
25-
GrB_Vector e,
26-
GrB_Index sink, // sink node
27-
GrB_Vector mask_vector,
28-
GrB_Vector active_set,
29-
GrB_Index *n_active,
30-
GrB_Index n,
24+
double *f, // total maxflow from src to sink
25+
GrB_Vector e, // excess vector, of type double
26+
GrB_Index sink, // sink node
27+
GrB_Vector src_and_sink, // mask vector, with just [src sink]
28+
// GrB_Vector active_set,
29+
GrB_Index *n_active, // # of active nodes
30+
// GrB_Index n,
3131
char *msg
3232
)
3333
{
34-
// f_T = e (T)
35-
double f_T = 0;
36-
GrB_Info info = GrB_Vector_extractElement(&f_T, e, sink); //if value at T
34+
// e_sink = e (sink)
35+
double e_sink = 0;
36+
GrB_Info info = GrB_Vector_extractElement(&e_sink, e, sink); //if value at sink
3737
GRB_TRY (info) ;
3838
if (info == GrB_SUCCESS)
3939
{
40-
// e(T) is present
41-
(*f) += f_T;
40+
// e(sink) is present
41+
(*f) += e_sink; /* FLOP */
4242
}
43-
GRB_TRY(GrB_select(active_set, mask_vector, NULL, GrB_VALUEGT_FP64, e, 0, GrB_DESC_RSC));
43+
44+
#if 0
45+
// e<![src,sink]> = select e where (e > 0)
46+
GRB_TRY(GrB_select(active_set, src_and_sink, NULL, GrB_VALUEGT_FP64, e, 0, GrB_DESC_RSC)); /* FLOP compare */
4447
// e = active_set : FIXME do we need this?
4548
GRB_TRY(GrB_assign(e, NULL, NULL, active_set, GrB_ALL, n, GrB_DESC_R));
4649
GRB_TRY(GrB_Vector_nvals(n_active, active_set));
50+
#else
51+
52+
// alternative:
53+
// e (src) = -1
54+
// e (sink) = -1
55+
// then remove src_and_sink
56+
57+
// TODO: what if e is tiny? Do we need a tol parameter,
58+
// and replace all "e > 0" and "r > 0" comparisons with
59+
// (e > tol), throughout the code?
60+
61+
// e<![src,sink]> = select e where (e > 0)
62+
GRB_TRY(GrB_select(e, src_and_sink, NULL, GrB_VALUEGT_FP64, e, 0, GrB_DESC_RSC)); /* FLOP compare */
63+
GRB_TRY(GrB_Vector_nvals(n_active, e));
64+
65+
#endif
4766
}
4867

4968
//------------------------------------------------------------------------------
@@ -61,14 +80,13 @@ static GrB_Info LG_augment_maxflow
6180
GrB_free(&theta); \
6281
GrB_free(&R); \
6382
GrB_free(&delta); \
64-
GrB_free(&delta); \
6583
GrB_free(&delta_vec); \
6684
GrB_free(&delta_mat); \
67-
GrB_free(&active_set); \
85+
/* GrB_free(&active_set); */ \
6886
GrB_free(&map); \
6987
GrB_free(&y); \
7088
GrB_free(&yd); \
71-
GrB_free(&mask_vector); \
89+
GrB_free(&src_and_sink); \
7290
GrB_free(&Jvec); \
7391
GrB_free(&GrB_Prune); \
7492
GrB_free(&GrB_UpdateFlows); \
@@ -90,16 +108,16 @@ static GrB_Info LG_augment_maxflow
90108
GrB_free(&GrB_InitBackwardFlows); \
91109
GrB_free(&GrB_CreateResidualForward); \
92110
GrB_free(&GrB_CreateResidualBackward); \
93-
GrB_free(&zero); \
111+
GrB_free(&zero); \
94112
GrB_free(&Re); \
95113
GrB_free(&invariant); \
96114
GrB_free(&GrB_InvariantCheck); \
97115
GrB_free(&check); \
98116
GrB_free(&GrB_extractYJ); \
99117
GrB_free(&extract_desc); \
100-
GrB_free(&residual_vec); \
101-
GrB_free(&GrB_MakeFlow); \
102-
GrB_free(&GrB_GetResidual); \
118+
GrB_free(&residual_vec); \
119+
GrB_free(&GrB_MakeFlow); \
120+
GrB_free(&GrB_GetResidual); \
103121
}
104122

105123

@@ -168,7 +186,7 @@ JIT_STR(void MF_CreateResidualBackward(MF_flowEdge *z, const double *y) {
168186
JIT_STR(void MF_RxdMult64(MF_resultTuple64 *z, const MF_flowEdge *x, GrB_Index ix,
169187
GrB_Index jx, const int64_t *y,
170188
GrB_Index iy, GrB_Index jy, const int64_t* theta) {
171-
double r = x->capacity - x->flow;
189+
double r = x->capacity - x->flow; /* FLOP */
172190
if(r > 0){
173191
z->d = *y;
174192
z->residual = r;
@@ -185,7 +203,7 @@ JIT_STR(void MF_RxdMult64(MF_resultTuple64 *z, const MF_flowEdge *x, GrB_Index i
185203
JIT_STR(void MF_RxdMult32(MF_resultTuple32 *z, const MF_flowEdge *x, GrB_Index ix,
186204
GrB_Index jx, const int32_t *y,
187205
GrB_Index iy, GrB_Index jy, const int32_t* theta) {
188-
double r = x->capacity - x->flow;
206+
double r = x->capacity - x->flow; /* FLOP */
189207
if(r > 0){
190208
z->d = *y;
191209
z->residual = r;
@@ -262,7 +280,7 @@ JIT_STR(void MF_extractFlow32(double *z, const MF_resultTuple32 *x)
262280
JIT_STR(void MF_updateFlow(MF_flowEdge *z,
263281
const MF_flowEdge *x, const double *y) {
264282
z->capacity = x->capacity;
265-
z->flow = x->flow + (*y);
283+
z->flow = x->flow + (*y); /* FLOP */
266284
}, GRB_UPDATEFLOWS_STR)
267285

268286

@@ -303,22 +321,22 @@ JIT_STR(void MF_extractYJ32(int32_t *z, const MF_resultTuple32 *x) {
303321

304322
JIT_STR(void MF_initForwardFlows(MF_flowEdge * z,
305323
const MF_flowEdge * x, const MF_flowEdge * y){
306-
z->flow = y->flow + x->flow;
324+
z->flow = y->flow + x->flow; /* FLOP */
307325
z->capacity = x->capacity;
308326
}, GRB_INITFLOWF_STR)
309327

310328

311329
JIT_STR(void MF_initBackwardFlows(MF_flowEdge * z,
312330
const MF_flowEdge * x, const MF_flowEdge * y){
313-
z->flow = x->flow - y->flow;
331+
z->flow = x->flow - y->flow; /* FLOP */
314332
z->capacity = x->capacity;
315333
}, GRB_INITFLOWB_STR)
316334

317335
JIT_STR(void MF_MxeMult64(MF_resultTuple64 * z, const MF_compareTuple64 * x,
318336
GrB_Index ix, GrB_Index jx,
319337
const double * y, GrB_Index iy,
320338
GrB_Index jy, const int64_t* theta){
321-
if(x->di == x->y_dmin && (*y) > 0){
339+
if(x->di == x->y_dmin && (*y) > 0){ /* FLOP compare */
322340
if(ix < jx){
323341
z->d = x->y_dmin;
324342
z->residual = x->residual;
@@ -330,7 +348,7 @@ JIT_STR(void MF_MxeMult64(MF_resultTuple64 * z, const MF_compareTuple64 * x,
330348
z->j = -1;
331349
}
332350
}
333-
else if(x->di == x->y_dmin - 1 && (*y) > 0){
351+
else if(x->di == x->y_dmin - 1 && (*y) > 0){ /* FLOP compare */
334352
z->d = INT64_MAX;
335353
z->residual = 0;
336354
z->j = -1;
@@ -351,7 +369,7 @@ JIT_STR(void MF_MxeMult32(MF_resultTuple32 * z, const MF_compareTuple32 * x,
351369
GrB_Index ix, GrB_Index jx,
352370
const double * y, GrB_Index iy,
353371
GrB_Index jy, const int32_t* theta){
354-
if(x->di == x->y_dmin && (*y) > 0){
372+
if(x->di == x->y_dmin && (*y) > 0){ /* FLOP compare */
355373
if(ix < jx){
356374
z->d = x->y_dmin;
357375
z->residual = x->residual;
@@ -363,7 +381,7 @@ JIT_STR(void MF_MxeMult32(MF_resultTuple32 * z, const MF_compareTuple32 * x,
363381
z->j = -1;
364382
}
365383
}
366-
else if(x->di == x->y_dmin - 1 && (*y) > 0){
384+
else if(x->di == x->y_dmin - 1 && (*y) > 0){ /* FLOP compare */
367385
z->d = INT32_MAX;
368386
z->residual = 0;
369387
z->j = -1;
@@ -526,7 +544,7 @@ JIT_STR(void MF_CheckInvariant32(bool *z, const int32_t *height,
526544

527545

528546
JIT_STR(void MF_getResidual(double * res, const MF_flowEdge * flow_edge){
529-
(*res) = flow_edge->capacity - flow_edge->flow;
547+
(*res) = flow_edge->capacity - flow_edge->flow; /* FLOP */
530548
}, GRB_GETRES_STR)
531549

532550

@@ -569,8 +587,8 @@ int LAGr_MaxFlow(double* f, LAGraph_Graph G, GrB_Index src, GrB_Index sink, char
569587
GrB_Vector d = NULL ;
570588

571589
//active_set and n_active
572-
GrB_Vector active_set = NULL ;
573-
GrB_Vector mask_vector = NULL ;
590+
//GrB_Vector active_set = NULL ;
591+
GrB_Vector src_and_sink = NULL ;
574592
GrB_Index n_active ;
575593

576594
//semiring and vectors for y<e, struct> = R x d
@@ -632,7 +650,7 @@ int LAGr_MaxFlow(double* f, LAGraph_Graph G, GrB_Index src, GrB_Index sink, char
632650
GRB_TRY(GrB_Matrix_nrows(&nrows, G->A));
633651
LG_ASSERT_MSG(nrows == n, GrB_INVALID_VALUE, "Matrix must be square");
634652
LG_ASSERT_MSG(src < n && src >= 0 && sink < n && sink >= 0,
635-
GrB_INVALID_VALUE, "S and T must be a value between [0, n)");
653+
GrB_INVALID_VALUE, "src and sink must be a value between [0, n)");
636654
LG_ASSERT_MSG(G->emin > 0, GrB_INVALID_VALUE,
637655
"the edge weights (capacities) must be greater than 0");
638656

@@ -692,15 +710,15 @@ int LAGr_MaxFlow(double* f, LAGraph_Graph G, GrB_Index src, GrB_Index sink, char
692710
GRB_TRY(GrB_assign(R, NULL, GrB_InitForwardFlows, Re, src, GrB_ALL, n, NULL));
693711
GRB_TRY(GrB_assign(R, NULL, GrB_InitBackwardFlows, Re, GrB_ALL, n, src, NULL));
694712

695-
//extract n_active from e masking T and S then assign to e
696-
GRB_TRY(GrB_Vector_new(&active_set, GrB_FP64, n));
697-
GRB_TRY(GrB_Vector_new(&mask_vector, GrB_BOOL, n)); //keep as bool?
698-
GRB_TRY (GrB_Vector_setElement (mask_vector, true, sink)) ;
699-
GRB_TRY (GrB_Vector_setElement (mask_vector, true, src)) ;
713+
//extract n_active from e masking sink and src then assign to e
714+
//GRB_TRY(GrB_Vector_new(&active_set, GrB_FP64, n));
715+
GRB_TRY(GrB_Vector_new(&src_and_sink, GrB_BOOL, n));
716+
GRB_TRY (GrB_Vector_setElement (src_and_sink, true, sink)) ;
717+
GRB_TRY (GrB_Vector_setElement (src_and_sink, true, src)) ;
700718

701-
// augment maxflow if the edge (S,T) exists
702-
LG_TRY (LG_augment_maxflow (f, e, sink, mask_vector,
703-
active_set, &n_active, n, msg)) ;
719+
// augment maxflow if the edge (src,sink) exists
720+
LG_TRY (LG_augment_maxflow (f, e, sink, src_and_sink,
721+
/* active_set, */ &n_active, /* n, */ msg)) ;
704722
//create flow vec
705723
GRB_TRY(GrB_Vector_new(&residual_vec, GrB_FP64, n));
706724

@@ -901,9 +919,9 @@ int LAGr_MaxFlow(double* f, LAGraph_Graph G, GrB_Index src, GrB_Index sink, char
901919
GRB_TRY(GrB_Matrix_new(&res_matT, GrB_FP64, n, n));
902920
GRB_TRY(GrB_Matrix_new(&res_mat, GrB_FP64, n, n));
903921
GRB_TRY(GrB_apply(res_mat, NULL,
904-
NULL, GrB_GetResidual, R, NULL)) ;
922+
NULL, GrB_GetResidual, R, NULL)) ; /* FLOP */
905923
GRB_TRY(GrB_select(res_mat, NULL,
906-
NULL, GrB_VALUEGT_FP64, res_mat, 0, NULL)) ;
924+
NULL, GrB_VALUEGT_FP64, res_mat, 0, NULL)) ; /* FLOP */
907925
GRB_TRY(GrB_transpose(res_matT, NULL, NULL, res_mat, NULL));
908926
LG_TRY(LAGraph_New(&res_graph,
909927
&res_matT, LAGraph_ADJACENCY_DIRECTED, msg));
@@ -913,15 +931,25 @@ int LAGr_MaxFlow(double* f, LAGraph_Graph G, GrB_Index src, GrB_Index sink, char
913931
LG_TRY(LAGr_BreadthFirstSearch(&lvl,
914932
NULL,
915933
res_graph, sink, msg));
916-
GRB_TRY(GrB_assign(d, mask_vector,
934+
935+
// FIXME: drop the src and sink
936+
// src_and_sink = [src sink]
937+
// d<![src,sink],struct> = lvl
938+
GRB_TRY(GrB_assign(d, src_and_sink,
917939
NULL, lvl, GrB_ALL, n, GrB_DESC_SC));
940+
// d<!lvl,struct> = n
918941
GRB_TRY(GrB_assign(d, lvl, NULL, n,
919942
GrB_ALL, n, GrB_DESC_SC));
943+
944+
// FIXME: simplify this to 1 call to GrB_assign [
945+
// e<!lvl,struct> = empty
920946
GRB_TRY(GrB_assign(e, lvl, NULL,
921947
-1, GrB_ALL, n, GrB_DESC_SC));
922948
GRB_TRY(GrB_select(e, NULL,
923-
NULL, GrB_VALUEGT_FP64,
949+
NULL, GrB_VALUEGT_FP64, /* FLOP */
924950
e, -1, NULL));
951+
// ]
952+
925953
GrB_free(&lvl);
926954
LG_TRY(LAGraph_Delete(&res_graph, msg));
927955
GRB_TRY(GrB_Vector_nvals(&n_active, e));
@@ -970,31 +998,34 @@ int LAGr_MaxFlow(double* f, LAGraph_Graph G, GrB_Index src, GrB_Index sink, char
970998

971999
//.min(flow_vec and e)
9721000
GRB_TRY(GrB_eWiseMult(delta_vec, NULL, NULL,
973-
GrB_MIN_FP64, residual_vec, e, NULL));
1001+
GrB_MIN_FP64, residual_vec, e, NULL)); /* FLOP */
9741002
GRB_TRY(GrB_apply(Jvec, NULL, NULL, GrB_extractYJ, y, NULL));
9751003
GRB_TRY(GrB_Matrix_clear(delta));
9761004
GRB_TRY(GxB_Matrix_build_Vector(delta, delta_vec,
9771005
Jvec, delta_vec, GxB_IGNORE_DUP, extract_desc));
9781006

979-
//make delta anti-symmetric
980-
GRB_TRY(GxB_eWiseUnion(delta_mat, NULL, NULL, GrB_MINUS_FP64,
1007+
//make delta anti-symmetric: delta_mat = (delta - delta')
1008+
GRB_TRY(GxB_eWiseUnion(delta_mat, NULL, NULL, GrB_MINUS_FP64, /* FLOP */
9811009
delta, zero, delta, zero, GrB_DESC_T1));
9821010

9831011
//update R
984-
GRB_TRY(GrB_eWiseMult(R, delta_mat, NULL, GrB_UpdateFlows,
1012+
// R<delta_mat> = UpdateFlows (R, delta_mat) using eWiseMult
1013+
GRB_TRY(GrB_eWiseMult(R, delta_mat, NULL, GrB_UpdateFlows, /* FLOP */
9851014
R, delta_mat, GrB_DESC_S));
9861015

9871016
//reduce delta_mat to delta_vec
1017+
// delta_vec = sum (delta_mat), summing up each row of delta_mat
9881018
GRB_TRY(GrB_reduce(delta_vec, NULL, NULL,
989-
GrB_PLUS_FP64, delta_mat, GrB_DESC_T0));
1019+
GrB_PLUS_FP64, delta_mat, GrB_DESC_T0)); /* FLOP */
9901020

9911021
//add to e
992-
GRB_TRY(GrB_assign(e, delta_vec, GrB_PLUS_FP64,
1022+
// e<delta_vec> += delta_vec
1023+
GRB_TRY(GrB_assign(e, delta_vec, GrB_PLUS_FP64, /* FLOP */
9931024
delta_vec, GrB_ALL, n, GrB_DESC_S));
9941025

9951026
// augment maxflow for all active nodes
996-
LG_TRY (LG_augment_maxflow (f, e, sink, mask_vector,
997-
active_set, &n_active, n, msg)) ;
1027+
LG_TRY (LG_augment_maxflow (f, e, sink, src_and_sink,
1028+
/* active_set, */ &n_active, /* n, */ msg)) ;
9981029

9991030
++iter;
10001031

0 commit comments

Comments
 (0)