@@ -149,6 +149,7 @@ GrB_Info LAGraph_CFL_reachability_multsrc
149149 GrB_Matrix MSrc ;
150150 GrB_Matrix M ;
151151 GrB_Matrix A ;
152+ GrB_Matrix B ;
152153 GrB_Vector a ;
153154 GrB_Index n ; // number of vertices in the graph
154155 GrB_Matrix identity_matrix = NULL ;
@@ -228,6 +229,7 @@ GrB_Info LAGraph_CFL_reachability_multsrc
228229
229230 GRB_TRY (GrB_Matrix_dup (& MSrc , TSrc [0 ]));
230231 GRB_TRY (GrB_Matrix_new (& A , GrB_BOOL , n , n ));
232+ GRB_TRY (GrB_Matrix_new (& B , GrB_BOOL , n , n ));
231233 GRB_TRY (GrB_Vector_new (& a , GrB_BOOL , n ));
232234 GRB_TRY (GrB_Matrix_new (& M , GrB_BOOL , n , n ));
233235
@@ -361,104 +363,25 @@ GrB_Info LAGraph_CFL_reachability_multsrc
361363 for (size_t i = 0 ; i < bin_rules_count ; i ++ ) {
362364 LAGraph_rule_WCNF bin_rule = rules [bin_rules [i ]];
363365
364- // #ifdef DEBUG_CFL_REACHABILITY
365- // printf("Rule: ");
366- // PRINT_RULE(bin_rule)
367- // #endif
368-
369- // #ifdef DEBUG_CFL_REACHABILITY
370- // printf("Before M = TSrc^A * T^B:\n");
371- // printf("TSrc^A:\n");
372- // PRINT_MATRIX(TSrc[bin_rule.nonterm]);
373- // printf("T^B:\n");
374- // PRINT_MATRIX(T[bin_rule.prod_A]);
375- // #endif
376-
377366 GRB_TRY (GrB_mxm (M , GrB_NULL , GrB_NULL , GxB_ANY_PAIR_BOOL ,
378- TSrc [bin_rule .nonterm ], T [bin_rule .prod_A ], GrB_NULL ));
379-
380- // #ifdef DEBUG_CFL_REACHABILITY
381- // printf("After M = TSrc^A * T^B:\n");
382- // printf("M:\n");
383- // PRINT_MATRIX(M);
384- // printf("Before T^A = T^A + M * T^C:\n");
385- // printf("T^A:\n");
386- // PRINT_MATRIX(T[bin_rule.nonterm]);
387- // printf("T^C:\n");
388- // PRINT_MATRIX(T[bin_rule.prod_B]);
389- // #endif
390-
391- GrB_BinaryOp acc_op = t_empty_flags [bin_rule .nonterm ] ? GrB_NULL : GxB_ANY_BOOL ;
392- GRB_TRY (GrB_mxm (T [bin_rule .nonterm ], GrB_NULL , acc_op , GxB_ANY_PAIR_BOOL ,
393- M , T [bin_rule .prod_B ], GrB_NULL ))
394-
395- // #ifdef DEBUG_CFL_REACHABILITY
396- // printf("After T^A = T^A + M * T^C:\n");
397- // printf("T^A:\n");
398- // PRINT_MATRIX(T[bin_rule.nonterm]);
399- // #endif
400-
401- // #ifdef DEBUG_CFL_REACHABILITY
402- // printf("Before TSrc^B = TSrc^B + TSrc^A:\n");
403- // printf("TSrc^A:\n");
404- // PRINT_MATRIX(TSrc[bin_rule.nonterm]);
405- // printf("TSrc^B:\n");
406- // PRINT_MATRIX(TSrc[bin_rule.prod_A]);
407- // #endif
367+ TSrc [bin_rule .nonterm ], T [bin_rule .prod_A ], GrB_NULL ));
408368
409- GRB_TRY (GrB_eWiseAdd (TSrc [bin_rule .prod_A ], GrB_NULL , GrB_NULL , GxB_ANY_BOOL ,
410- TSrc [bin_rule .prod_A ], TSrc [bin_rule .nonterm ], GrB_NULL ));
411369
412- // #ifdef DEBUG_CFL_REACHABILITY
413- // printf("After TSrc^B = TSrc^B + TSrc^A:\n");
414- // printf("TSrc^B:\n");
415- // PRINT_MATRIX(TSrc[bin_rule.prod_A]);
416- // #endif
370+ GRB_TRY (GrB_mxm (B , GrB_NULL , GrB_NULL , GxB_ANY_PAIR_BOOL , M , T [bin_rule .prod_B ], GrB_NULL ));
371+
372+ GRB_TRY (GrB_eWiseAdd (T [bin_rule .nonterm ], GrB_NULL , GrB_NULL , GxB_ANY_BOOL , T [bin_rule .nonterm ], B , GrB_NULL ));
417373
418- // Update source vertices matrix to find appropriate paths only
419374
420- // #ifdef DEBUG_CFL_REACHABILITY
421- // printf("Before A = dest(M)\n");
422- // printf("M:\n");
423- // PRINT_MATRIX(M);
424- // #endif
375+ GRB_TRY (GrB_eWiseAdd (TSrc [bin_rule .prod_A ], GrB_NULL , GrB_NULL , GxB_ANY_BOOL ,
376+ TSrc [bin_rule .prod_A ], TSrc [bin_rule .nonterm ], GrB_NULL ));
425377
378+ // Update source vertices matrix to find appropriate paths only
426379 // M[i, j] == 1 => A[j, j] == 1
427380 GRB_TRY (GrB_vxm (a , GrB_NULL , GrB_NULL , GxB_ANY_PAIR_BOOL , ones_vec , M , GrB_NULL ));
428381 GrB_Matrix_free (& A );
429- GRB_TRY (GrB_Matrix_diag (& A , a , 0 ));
430-
431- // #ifdef DEBUG_CFL_REACHABILITY
432- // printf("After A = dest(M)\n");
433- // printf("a:\n");
434- // PRINT_VECTOR(a);
435- // printf("A:\n");
436- // PRINT_MATRIX(A);
437- // #endif
438-
439- // #ifdef DEBUG_CFL_REACHABILITY
440- // printf("Before TSrc^C = TSrc^c + A\n");
441- // printf("TSrc^C:\n");
442- // PRINT_MATRIX(TSrc[bin_rule.prod_B])
443- // printf("A:\n");
444- // PRINT_MATRIX(A)
445- // #endif
446-
382+ GRB_TRY (GrB_Matrix_diag (& A , a , 0 );
447383 GRB_TRY (GrB_eWiseAdd (TSrc [bin_rule .prod_B ], GrB_NULL , GrB_NULL , GxB_ANY_BOOL ,
448- TSrc [bin_rule .prod_B ], A , GrB_NULL ));
449-
450- // #ifdef DEBUG_CFL_REACHABILITY
451- // printf("After TSrc^C = TSrc^c + A\n");
452- // printf("TSrc^C:\n");
453- // PRINT_MATRIX(TSrc[bin_rule.prod_B])
454- // #endif
455-
456- // #ifdef DEBUG_CFL_REACHABILITY
457- // GxB_Matrix_iso(&iso_flag, T[bin_rule.nonterm]);
458- // printf("[TERM1 TERM2] MULTIPLY, S: %d, A: %d, B: %d, "
459- // "I: %ld (ISO: %d)\n",
460- // bin_rule.nonterm, bin_rule.prod_A, bin_rule.prod_B, i, iso_flag);
461- // #endif
384+ TSrc [bin_rule .prod_B ], A , GrB_NULL )));
462385
463386 // Check if any of the matrices changed. If not, job is done.
464387 GrB_Index nnz_T , nnz_TSrc_B , nnz_TSrc_C ;
@@ -501,5 +424,6 @@ GrB_Info LAGraph_CFL_reachability_multsrc
501424 GRB_TRY (GrB_Matrix_dup (output , MSrc ));
502425
503426 LG_FREE_ALL ;
427+
504428 return GrB_SUCCESS ;
505429}
0 commit comments