@z 
This file was created by Aleksandar Donev as part of the Network Optimization project. 
Feel free to use any portion of it and contact me at donev@pa.msu.edu 
@x 
\Title{Graph and Network Algorithms from Public FORTRAN 77 codes}
\author{Aleksandar Donev}
\date{July 2001}
\maketitle

@*0 Module |Graph_Algorithms_F77|.

This file contains interfaces to a shortes-path Dijkstra public-domain F77 code, as well as an interface to a min-cut F77 code. These codes are in the file |"GraphsF77.f90"|, where there is a header giving the original authors. These both calculate a given number of shortest paths or a given number of recursive cuts, to be used in comparisons with SSCNO solutions.

This documentation is not finished!

@
@a
MODULE Graph_Algorithms_F77 @;
   USE Precision @;
   USE Error_Handling @;
   USE System_Monitors @;
   USE Sorting_Ranking @;
   USE Graph_Algorithms @;   
   USE Network_Matrix_Operations @;
   IMPLICIT NONE @; 
   PUBLIC :: DijkstraShortestPaths, LabelingMinCut @;
 
   PRIVATE @;

	INTERFACE @; // Imported codes	
	   SUBROUTINE ShortestPaths_F77(A, ND, LNGT, D, P, ARR1, ARR2, NMAX, MMAX, N, R, ORD) @;  
	      INTEGER A, ND, LNGT, D, P, ARR1, ARR2, NMAX, MMAX, N, R @;
	      DIMENSION A (:), D (:), P (:), ARR1 (:), ARR2 (:), ND (:), LNGT (:) @;
	      CHARACTER ORD @;
	   END SUBROUTINE ShortestPaths_F77 @;
      SUBROUTINE MinCut_F77(DNFROM, DNTO, DNGCAP, NMASK, NODES, NARCS, &
                      ISOURC, ISINK, MAXFLO, NNCUT, NACUT, IRETN) @;
         INTEGER  :: DNFROM (:), DNTO (:), DNGCAP(:) @;
         INTEGER  :: NODES, NARCS, ISOURC, ISINK, MAXFLO, NNCUT, NACUT, IRETN @;
         LOGICAL(1) :: NMASK(:) @;
	   END SUBROUTINE MinCut_F77 @;
	END INTERFACE @;            
                                                 
CONTAINS @;

   @<DijkstraShortestPaths@>@;
   @<LabelingMinCut@>

END MODULE Graph_Algorithms_F77 @;   

@*1 Shortest Paths.

Calculates the Dijkstra shortest tree from S to all neighbours of T and selects and marks the arcs on the first |n_paths| shortest paths:

@
@<DijkstraShortestPaths@>=
SUBROUTINE DijkstraShortestPaths(node_offset, heads_tails, special_heads_tails, arcs_lengths, &
                                 source_sink, shortest_paths, n_paths, shortest_lengths, queue_type) @;

   @%% USE Network_Data_Structures, ONLY : nodes_coords @;
   @%% USE Network_Graphics @;                                  
   IMPLICIT NONE @;
   INTEGER(KIND=i_wp), INTENT(IN) :: node_offset @;                               
   INTEGER(KIND=i_wp), DIMENSION(:,:), INTENT(IN) :: heads_tails, special_heads_tails @;
   REAL(KIND=r_wp), DIMENSION(:), INTENT(IN) :: arcs_lengths @; // Weights on regular arcs    
   INTEGER(KIND=i_wp), DIMENSION(2), INTENT(IN) :: source_sink @; // S and T
   INTEGER(KIND=i_byte), DIMENSION(-node_offset:), INTENT(OUT) :: shortest_paths @; // The shortest path a node belongs to
   INTEGER(KIND=i_byte), INTENT(IN) :: n_paths @; // The number of first shortest paths to output (up to 127)
   REAL(KIND=r_wp), DIMENSION(:), INTENT(OUT), OPTIONAL :: shortest_lengths @; // Ratio of shortest path length to this path's length
   CHARACTER, INTENT(IN), OPTIONAL :: queue_type @; // 'Unordered', 'Heap' or 'List'

   INTEGER(KIND=i_wp), DIMENSION(:), ALLOCATABLE :: &
                         neighbours, my_neighbours, incident_arcs, &
                         parents, distances, temp_array_1, temp_array_2, edge_lengths, &
                         paths_lengths, paths_ordering @;
   INTEGER(KIND=i_wp) :: n_nodes, n_special_nodes, n_arcs, n_special_arcs, &
                         n_edges, n_vertices @; // Counters
   INTEGER(KIND=i_wp) :: node, arc, head, tail, special_node, special_arc, counter, &
                         min_distance, max_distance, vertex, edge, &
                         max_edge_weight, root, s_or_t, ST_degrees[2], &
                         source_vertex, sink_vertex @;
   INTEGER(KIND=i_byte) :: path @;                         
   REAL(KIND=r_wp) :: max_weight @;                                                   
   INTEGER :: alloc_status, status @;      
   CHARACTER :: queue @; 

   queue='H' @; // Use a heap by default
   IF(PRESENT(queue_type)) @~ queue=queue_type @;        
   
   n_special_nodes=-_LBOUND(shortest_paths, i_wp, DIM=1) @;
   n_nodes=_UBOUND(shortest_paths, i_wp, DIM=1) @; // The number of nodes, $n$ @;
   n_special_arcs=_SIZE(special_heads_tails, i_wp, DIM=2) @;
   n_arcs=_SIZE(heads_tails, i_wp, DIM=2) @; // $m$

   IF(ANY(source_sink>=0)) THEN @;
      // Here the ST nodes must be special
      CALL NonCriticalError("Source or sink is not a special node", "DijkstraShortestPaths") @;
      RETURN @;
   END IF @;
   n_vertices=n_nodes+2 @; // For now I add only the ST special nodes
   root=n_nodes+1 @; // Root of spanning tree
   n_edges=2*(n_arcs+n_special_arcs) @; // Maximum possible number of edges

   _AllocateArray(my_neighbours, 0, n_vertices+1, 1_i_wp, .NOT.ALLOCATED, "DijkstraShortestPaths") @;
   _AllocateArray(distances, 0, n_vertices+1, 1_i_wp, .NOT.ALLOCATED, "DijkstraShortestPaths") @;
   _AllocateArray(parents, 0, n_vertices+1, 1_i_wp, .NOT.ALLOCATED, "DijkstraShortestPaths") @;

   _AllocateArray(neighbours, 0, n_edges, 1_i_wp, .NOT.ALLOCATED, "DijkstraShortestPaths") @;               
   _AllocateArray(incident_arcs, 0, n_edges, 1_i_wp, .NOT.ALLOCATED, "DijkstraShortestPaths") @;
   _AllocateArray(edge_lengths, 0, n_edges, 1_i_wp, .NOT.ALLOCATED, "DijkstraShortestPaths") @;        

   @%% WRITE(*,*) n_nodes, n_special_nodes, n_arcs, n_special_arcs, n_vertices, n_edges
   CALL CreateAdjacencyArrays(arc_offset=-1, node_offset=-1, &
               heads_tails=heads_tails, my_neighbours=my_neighbours[0:n_nodes+1], &
               neighbours=neighbours, incident_arcs=incident_arcs, n_edges=n_edges) @;
       
   max_weight=MAXVAL(arcs_lengths) @; // Maybe logarithmic scaling?:
   max_edge_weight=HUGE(1_i_wp)/(10*n_vertices) @;
   DO edge=1, n_edges @; // Normalization of edge weights
      arc=incident_arcs[edge] @;
      edge_lengths[edge]=INT(@E REAL(max_edge_weight)*arcs_lengths[arc]/max_weight, i_wp)+1_i_wp @;      
   END DO @;
    
   ST_degrees=0 @; // The degrees of S and T
   STArcs: DO s_or_t=1,2 @;
      my_neighbours[n_nodes+s_or_t]=n_edges+1 @;
	   DO special_arc=1, n_special_arcs @; // We add the ST arcs now
	      IF((special_heads_tails[1,special_arc]==source_sink[s_or_t]).AND. &
	         (special_heads_tails[2,special_arc]>0)) THEN @;
	         ST_degrees[s_or_t]=ST_degrees[s_or_t]+1 @;
	         n_edges=n_edges+1 @;
            neighbours[n_edges]=special_heads_tails[2,special_arc] @;
            edge_lengths[n_edges]=1 @;
            incident_arcs[n_edges]=-special_arc @;
         END IF @;   
         IF((special_heads_tails[2,special_arc]==source_sink[s_or_t]).AND. &
            (special_heads_tails[1,special_arc]>0)) THEN @;       
	         ST_degrees[s_or_t]=ST_degrees[s_or_t]+1 @;
	         n_edges=n_edges+1 @;
            neighbours[n_edges]=special_heads_tails[1,special_arc] @;
            edge_lengths[n_edges]=1 @;  
            incident_arcs[n_edges]=-special_arc @;          
         END IF @;   
      END DO @;      
   END DO STArcs @;
   my_neighbours[n_vertices+1]=n_edges+1 @; // Sentinel value
   
   @#if 0               
	   WRITE(*,"(25I5)") (node, node=1, n_edges) @;
	   WRITE(*,"(25I5)") my_neighbours[1:n_vertices] @;   
	   WRITE(*,"(25I5)") neighbours[1:n_edges] @;    
	   WRITE(*,"(25I5)") incident_arcs[1:n_edges] @;
	   WRITE(*,"(25I5)") edge_lengths[1:n_edges] @; 
	   PAUSE @;  
   @#endif
     
   _AllocateArray(temp_array_1, 0, n_vertices+1, 1_i_wp, .NOT.ALLOCATED, "DijkstraShortestPaths") @;
   _AllocateArray(temp_array_2, 0, n_vertices+1, 1_i_wp, .NOT.ALLOCATED, "DijkstraShortestPaths") @;
   
   CALL ShortestPaths_F77 (A=my_neighbours[1:n_vertices], ND=neighbours[1:n_edges], &
         LNGT=edge_lengths[1:n_edges], D=distances[1:n_vertices], &
         P=parents[1:n_vertices], ARR1=temp_array_1[1:n_vertices], &
         ARR2=temp_array_2[1:n_vertices], NMAX=n_vertices, MMAX=n_edges, &
         N=n_nodes+1, R=root, ORD=queue) @;
         
   _DeallocateArray(temp_array_1, 1_i_wp, ALLOCATED,"DijkstraShortestPaths") @;
   _DeallocateArray(temp_array_2, 1_i_wp, ALLOCATED,"DijkstraShortestPaths") @;               

   _AllocateArray(paths_lengths, 1, ST_degrees[2], 1_i_wp, .NOT.ALLOCATED, "DijkstraShortestPaths") @;
   _AllocateArray(paths_ordering, 1, ST_degrees[2], 1_i_wp, .NOT.ALLOCATED, "DijkstraShortestPaths") @;

   source_vertex=n_nodes+1 @;
   sink_vertex=n_nodes+2 @;
   paths_lengths=distances(neighbours(my_neighbours[sink_vertex]:(my_neighbours[sink_vertex+1]-1))) @;
   CALL QuickRank(array=paths_lengths, permutation=paths_ordering) @;
   distances[sink_vertex]=paths_lengths(paths_ordering[1]) @;

   @#if 0
	   CALL InitNetworkGraphics(file="D",file_type="CONS",  &                       
                       page_size=(/5000,5000/), label_format="5N0", &
                       color_table="RAIN", colorbar_position="Horizontal", &
                       axis_labels_format=(/"NONE","NONE","NONE","NONE"/) ) @;                                                                      
	   CALL PlotNetwork2D(heads_tails=heads_tails[:,1:n_arcs], &
	                     node_offset=-1, node_coords=nodes_coords[:,1:n_nodes], &
	                     node_values=@E REAL(distances[1:n_nodes], r_wp), &
	                     arc_values=arcs_lengths[1:n_arcs], &	                       
	                     node_size_range=(/-HUGE(1.0_r_wp)/100, HUGE(1.0_r_wp)/20/), &
	                     vector_type=0) @;	                     
      CALL EndNetworkGraphics() @;	                     
    
	   WRITE(*,*) "___ " @;
	   @%% WRITE(*,*) "All distances:", distances[1:n_vertices] @;
	   WRITE(*,*) "T nodes:", neighbours(my_neighbours[sink_vertex]:(my_neighbours[sink_vertex+1]-1)) @;
	   WRITE(*,*) "Distances of T nodes:", paths_lengths @;
	   WRITE(*,*) "Ordered distances of T nodes:", paths_lengths(paths_ordering) @;
	   WRITE(*,*) "Order of T nodes:", neighbours(my_neighbours[sink_vertex]+paths_ordering-1) @;
	   WRITE(*,*) "___ " @;
	   PAUSE @;
   @#endif

   min_distance=paths_lengths(paths_ordering[1]) @;
   max_distance=paths_lengths(paths_ordering[ST_degrees(2)]) @;
   WRITE(message_log_unit,*) "The range of distances of T nodes is:", min_distance, max_distance @;
      
   shortest_paths=0 @; // Sentinel value
   SelectedPaths: DO path=1_i_byte, INT(MIN(INT(n_paths, i_wp), ST_degrees[2]), i_wp)  @;      
      vertex=neighbours(my_neighbours[sink_vertex]+paths_ordering[path]-1) @;
      IF(PRESENT(shortest_lengths)) @~ shortest_lengths[path]=@E REAL(min_distance)/@E REAL(distances[vertex]) @;
      MarkPath: DO @;
         counter=counter+1 @;
         node=vertex @; // The numbering is the same         
         IF((shortest_paths[node]!=0).OR.(vertex==source_vertex)) @~ EXIT MarkPath @;
         shortest_paths[node]=path @; // Mark this node on the path
         vertex=parents[vertex] @; // Go to the next vertex
      END DO MarkPath @;
   END DO SelectedPaths @;        
   @%% WRITE(*,*) "Node labels:", shortest_paths @;
      
   _DeallocateArray(paths_lengths, 1_i_wp, ALLOCATED,"DijkstraShortestPaths") @;
   _DeallocateArray(paths_ordering, 1_i_wp, ALLOCATED,"DijkstraShortestPaths") @;
   
   _DeallocateArray(neighbours, 1_i_wp, ALLOCATED,"DijkstraShortestPaths") @;
   _DeallocateArray(parents, 1_i_wp, ALLOCATED,"DijkstraShortestPaths") @;
   _DeallocateArray(distances, 1_i_wp, ALLOCATED,"DijkstraShortestPaths") @;
               
   _DeallocateArray(my_neighbours, 1_i_wp, ALLOCATED,"DijkstraShortestPaths") @;
   _DeallocateArray(incident_arcs, 1_i_wp, ALLOCATED,"DijkstraShortestPaths") @;
   _DeallocateArray(edge_lengths, 1_i_wp, ALLOCATED,"DijkstraShortestPaths") @;   
         
END SUBROUTINE DijkstraShortestPaths @;

@*1 Min-Cut/Max-Flow Algorithm.

Calculates |n_cuts| recursive minimal cuts, starting with an ST cut and then working on the resulting halves.

@
@<LabelingMinCut@>=
SUBROUTINE LabelingMinCut(node_offset, arc_offset, heads_tails, arcs_capacities, &
                          source_sink, nodes_cuts, arcs_cuts, n_cuts, cuts_capacities) @;

   @%% USE Network_Data_Structures, ONLY : nodes_coords @;
   @%% USE Network_Graphics @;                                  
   IMPLICIT NONE @;
   INTEGER(KIND=i_wp), INTENT(IN) :: arc_offset, node_offset @;                               
   INTEGER(KIND=i_wp), DIMENSION(:,-arc_offset:), INTENT(IN) :: heads_tails @;
   REAL(KIND=r_wp), DIMENSION(1:), INTENT(IN) :: arcs_capacities @; // Capacities of regular arcs {\em only}    
   INTEGER(KIND=i_wp), DIMENSION(2), INTENT(IN) :: source_sink @; // S and T
   INTEGER(KIND=i_byte), DIMENSION(-node_offset:), INTENT(OUT) :: nodes_cuts @; // The cut tree node
   INTEGER(KIND=i_byte), DIMENSION(:), INTENT(OUT) :: arcs_cuts @;
   INTEGER(KIND=i_byte), INTENT(IN) :: n_cuts @; // The number of recursive cut calculations (less then 7: |2^n_cuts<HUGE(1_i_wp)|)
   REAL(KIND=r_wp), DIMENSION(:), INTENT(OUT), OPTIONAL :: cuts_capacities @; // Ratio of largest cut to this cut's capacity

   INTEGER(KIND=i_byte), DIMENSION(:), ALLOCATABLE :: cuts_renumbering @;
   INTEGER(KIND=i_wp), DIMENSION(:), ALLOCATABLE :: cuts_reordering @;   
   INTEGER(KIND=i_wp) :: n_nodes, n_special_nodes, n_arcs, n_special_arcs, &
                         sub_n_nodes, sub_n_arcs, n_vertices, n_edges @; // Counters
   INTEGER(KIND=i_wp) :: node, arc, head, tail, special_node, special_arc, counter, &
                         min_distance, max_distance, vertex, edge, &
                         max_edge_weight, root, s_or_t, ST_degrees[2], &
                         source_vertex, sink_vertex @;
   INTEGER(KIND=i_byte) :: cut, low_cut, high_cut, parent_cut, cut_level, parent_label @;                         
   REAL(KIND=r_wp) :: max_capacity, normalization @;                                                   
   INTEGER :: alloc_status, status @;   
   @%% trace_allocations=.TRUE. @;
   
   n_special_nodes=-_LBOUND(nodes_cuts, i_wp, DIM=1) @;
   n_nodes=_UBOUND(nodes_cuts, i_wp, DIM=1) @; // The number of nodes, $n$ @;
   n_special_arcs=-_LBOUND(heads_tails, i_wp, DIM=2) @;
   n_arcs=_UBOUND(heads_tails, i_wp, DIM=2) @;    
   @%% WRITE(*,*) "Top: ", n_nodes, n_special_nodes, n_arcs, n_special_arcs @;

   IF(ANY(source_sink>=0)) THEN @;
      // Here the ST nodes must be special
      CALL NonCriticalError("Source or sink is not a special node", "LabelingMinCut") @;
      RETURN @;
   END IF @;

   max_capacity=MAXVAL(arcs_capacities) @; // Maybe logarithmic scaling?:
   normalization= @E REAL(HUGE(1_i_wp)/10/n_nodes, r_wp)/max_capacity @;
   
   nodes_cuts[1:n_nodes]=0 @;
   nodes_cuts[-n_special_nodes:0]=-2 @; // A dummy value
   nodes_cuts[source_sink(1)]=-1 @; // The label for S
   nodes_cuts[source_sink(2)]=HUGE(1_i_byte) @; // The label for T
   arcs_cuts=0 @;
   
   cut_level=0 @;
   low_cut=0 @; // Initalize
   high_cut=0 @; // Initially
   cut=0 @;
   cuts_capacities=0.0_r_wp @;
   RecursiveCuts: DO @;
      cut_level=cut_level+1 @;
      nodes_cuts[1:n_nodes]=2_i_byte*nodes_cuts[1:n_nodes]+1_i_byte @;

      DO parent_cut=low_cut, high_cut @;

         cut=cut+1 @;
         IF(cut>n_cuts) @~ EXIT RecursiveCuts @;
         parent_label=INT(2*parent_cut+1,i_byte) @; 
         @%% WRITE(*,*) "Working on parent cut:", parent_cut, " labeled", parent_label, &
         @%%   " label range:", MINVAL(nodes_cuts[1:n_nodes]), MAXVAL(nodes_cuts[1:n_nodes]) @;
  
         CALL SubMinCut(node_offset=n_special_nodes, arc_offset=n_special_arcs, heads_tails=heads_tails, &
                  arcs_capacities=arcs_capacities, capacity_normalization=normalization, &
                  parent_cut=parent_label, cuts_labels=nodes_cuts, cut_capacity=cuts_capacities[cut]) @;   

	      MarkArcs: DO arc=1, n_arcs @; // Mark the arcs in this cut
	         head=heads_tails[1,arc] @;
	         tail=heads_tails[2,arc] @;            
	         IF(((nodes_cuts[head]==parent_label).AND.(nodes_cuts[tail]==parent_label+1)).OR. &
	            ((nodes_cuts[head]==parent_label+1).AND.(nodes_cuts[tail]==parent_label))) THEN @;
	            arcs_cuts[arc]=cut @;
	            @%% WRITE(*,*) "Arc, cut, capacity:", arc, cut, cuts_capacities[cut] @;
            END IF @;     
	      END DO MarkArcs @;
	      @%% WRITE(*,*) "old cut, capacity:", cut, cuts_capacities[cut] @;
                             
      END DO @;

      low_cut=high_cut+1 @;
      high_cut=low_cut+2^cut_level-1 @;
   END DO RecursiveCuts @;   
   nodes_cuts[-n_special_nodes:0]=0 @; // Return dummy values
   
   // Now order the cuts in order of capacity:
  _AllocateArray(cuts_reordering, 1, n_cuts, 1_i_wp, .NOT.ALLOCATED, "LabelingMinCut") @;
  _AllocateArray(cuts_renumbering, 1, n_cuts, 1_i_byte, .NOT.ALLOCATED, "LabelingMinCut") @;  
   CALL QuickRank(array=cuts_capacities[1:n_cuts], permutation=cuts_reordering) @;
   
   DO cut=1, n_cuts @; // Reorder the cuts sizes from largest to smallest
      cuts_renumbering(cuts_reordering[cut])=n_cuts-cut+1 @;
   END DO @;   
   cuts_capacities[1:n_cuts]=cuts_capacities(cuts_reordering[n_cuts:1:(-1)])/cuts_capacities(cuts_reordering[n_cuts]) @;      

   @#if 0
	   DO cut=1, n_cuts @; // Reorder the cuts sizes from largest to smallest
	      WRITE(*,*) "old cut, new cut, capacity:", cut, cuts_renumbering(cut), cuts_capacities(cuts_renumbering[cut]) @; 
	   END DO @;   
   @#endif
   
	DO arc=1, n_arcs @; // Mark the arcs in this cut
      IF(arcs_cuts[arc]>0) THEN @;
         @%% WRITE(*,*) "arc, old cut, new cut, capacity", arc, arcs_cuts[arc], &
         @%%             cuts_renumbering(arcs_cuts[arc]), cuts_capacities(cuts_renumbering(arcs_cuts[arc])) @;      
         arcs_cuts[arc]=cuts_renumbering(arcs_cuts[arc]) @; 
      END IF @;   
	END DO @;
	
   _DeallocateArray(cuts_reordering, 1_i_wp, ALLOCATED, "LabelingMinCut") @;
   _DeallocateArray(cuts_renumbering, 1_i_byte, ALLOCATED, "LabelingMinCut") @;     
   @%% trace_allocations=.FALSE. @;
   @%% STOP
                                 
CONTAINS @;

   SUBROUTINE SubMinCut(node_offset, arc_offset, heads_tails, arcs_capacities, capacity_normalization, &
                        parent_cut, cuts_labels, cut_capacity) @;
	   IMPLICIT NONE @;                  
      INTEGER(KIND=i_wp), INTENT(IN) :: node_offset, arc_offset @;      	    
	   INTEGER(KIND=i_wp), DIMENSION(:,-arc_offset:), INTENT(IN) :: heads_tails @;  
	   REAL(KIND=r_wp), DIMENSION(1:), INTENT(IN) :: arcs_capacities @; // Arc capacities of {\regular arcs} only
	   REAL(KIND=r_wp), INTENT(IN), OPTIONAL :: capacity_normalization @;   
	   INTEGER(KIND=i_byte), INTENT(IN) :: parent_cut @; // S nodes will get label |parent_cut|, T nodes |parent_cut+1|
	   INTEGER(KIND=i_byte), DIMENSION(-node_offset:), INTENT(INOUT) :: cuts_labels @; // The labels of nodes
	   @%% INTEGER(KIND=i_byte), DIMENSION(:), INTENT(INOUT) :: arcs_cuts @;
      REAL(KIND=r_wp), INTENT(OUT) :: cut_capacity @;	   

      INTEGER(KIND=i_wp), DIMENSION(:), ALLOCATABLE :: heads, tails, capacities @;	   
      INTEGER(KIND=i_wp), DIMENSION(:), ALLOCATABLE :: nodes_renumbering, nodes_reordering @;
      LOGICAL(KIND=l_short), DIMENSION(:), ALLOCATABLE :: ST_mask @;
   
      INTEGER(KIND=i_wp) :: n_nodes, n_S_nodes, n_T_nodes, n_arcs, n_special_arcs, n_S_arcs, n_T_arcs, &
                            n_cut_arcs, n_vertices, n_edges, source_vertex, sink_vertex @;      
      INTEGER(KIND=i_wp) :: node, arc, edge, vertex, head, tail, counter, i_temp @;
      INTEGER(KIND=i_wp), PARAMETER :: inf_capacity=HUGE(1_i_wp)/10 @; // This will cause cuts 
      INTEGER(KIND=i_byte) :: cut @;
      REAL(KIND=r_wp) :: normalization @;

      normalization=1.0_r_wp @;
      IF(PRESENT(capacity_normalization)) @~ normalization=capacity_normalization @;
      
      n_special_nodes=-_LBOUND(cuts_labels, i_wp, DIM=1) @;
      n_nodes=_UBOUND(cuts_labels, i_wp, DIM=1) @;    
	   n_special_arcs=-_LBOUND(heads_tails, i_wp, DIM=2) @;
	   n_arcs=_UBOUND(heads_tails, i_wp, DIM=2) @; 
      @%% WRITE(*,*) "Inside: ", n_nodes, n_special_nodes, n_arcs, n_special_arcs @;

      _AllocateNodalArray(nodes_renumbering, 1_i_wp, .NOT.ALLOCATED, "LabelingMinCut") @;  

      nodes_renumbering=0 @; // For masking purposes
      counter=0 @;
      DO node=-n_special_nodes, n_nodes @;
         IF(cuts_labels[node]==parent_cut) THEN @;
            counter=counter+1 @;
            nodes_renumbering[node]=counter @; // The new node number
         END IF @;
      END DO @;
      n_vertices=counter+2 @;
      source_vertex=n_vertices-1 @;
      sink_vertex=n_vertices @;

      n_edges=0 @;
      DO arc=-n_special_arcs, n_arcs @; // Count how many arcs enter in this subgraph
         head=heads_tails[1,arc] @;
         tail=heads_tails[2,arc] @;            
         // Each arc can contribute to at most two edges
         IF((cuts_labels[head]==parent_cut).AND.(cuts_labels[tail]==parent_cut)) THEN @;
            n_edges=n_edges+2 @;
         ELSE IF((cuts_labels[head]==parent_cut).OR.(cuts_labels[tail]==parent_cut)) THEN @;
            n_edges=n_edges+1 @;
         END IF @;      
      END DO @;
      
      _AllocateArray(heads, 1, n_edges, 1_i_wp, .NOT.ALLOCATED, "LabelingMinCut") @;
      _AllocateArray(tails, 1, n_edges, 1_i_wp, .NOT.ALLOCATED, "LabelingMinCut") @;   
      _AllocateArray(capacities, 1, n_edges, 1_i_wp, .NOT.ALLOCATED, "LabelingMinCut") @;  
      _AllocateArray(ST_mask, 1, n_vertices, .TRUE._l_short, .NOT.ALLOCATED, "LabelingMinCut") @;          

      @#if 1
      edge=0 @;
      n_S_arcs=0 @;
      n_T_arcs=0 @;
      DO arc=-n_special_arcs, n_arcs @;
         head=heads_tails[1,arc] @;
         tail=heads_tails[2,arc] @;            

         IF((cuts_labels[head]==parent_cut).AND.(cuts_labels[tail]==parent_cut)) THEN @;

            edge=edge+1 @;
            heads[edge]=nodes_renumbering[head] @;
            tails[edge]=nodes_renumbering[tail] @;            
            capacities[edge]=INT(normalization*arcs_capacities[arc],i_wp) @;
            edge=edge+1 @; // We need to add this arc twice as an edge
            heads[edge]=nodes_renumbering[tail] @;
            tails[edge]=nodes_renumbering[head] @;
            IF(arc>0) THEN @;
               capacities[edge]=INT(normalization*arcs_capacities[arc],i_wp) @;
            ELSE @;
               capacities[edge]=inf_capacity @; // Infinite capacity for special arcs
            END IF @;   

         ELSE IF(cuts_labels[head]==parent_cut) THEN @;

            @%% WRITE(*,*) "Adding ST edge:", edge, head, tail @;
            edge=edge+1 @;
            IF(cuts_labels[tail]<parent_cut) THEN @; // Connect this to S
               n_S_arcs=n_S_arcs+1 @;
               tails[edge]=source_vertex @; // This arc emanates from S
               heads[edge]=nodes_renumbering[head] @; // And terminates at the node                       
            ELSE @; // Connect it to T
               n_T_arcs=n_T_arcs+1 @;
               heads[edge]=sink_vertex @;
               tails[edge]=nodes_renumbering[head] @;
            END IF @;
            capacities[edge]=inf_capacity @; // Infinite capacity

         ELSE IF(cuts_labels[tail]==parent_cut) THEN @;
         
            @%% WRITE(*,*) "Adding ST edge:", edge, head, tail @;
            edge=edge+1 @;
            IF(cuts_labels[head]<parent_cut) THEN @; // Connect this to S
               n_S_arcs=n_S_arcs+1 @;
               tails[edge]=source_vertex @; // This arc emanates from S
               heads[edge]=nodes_renumbering[tail] @; // And terminates at the node                       
            ELSE @; // Connect it to T
               n_T_arcs=n_T_arcs+1 @;
               heads[edge]=sink_vertex @;
               tails[edge]=nodes_renumbering[tail] @;
            END IF @;
            capacities[edge]=inf_capacity @; // Infinite capacity

         END IF @;
         
      END DO @;
      @%% WRITE(*,*) "For parent ", parent_cut, " n_ST_arcs=", n_S_arcs, n_T_arcs @;
      @#endif
         
      @#if 0
         WRITE(*,*) "ST:", source_vertex, sink_vertex @;         
         WRITE(*,*) tails @;
         WRITE(*,*) heads @;
         WRITE(*,*) capacities @;
         @%% PAUSE @;
      @#endif

      @#if 1
      CALL MinCut_F77(DNFROM=tails, DNTO=heads, DNGCAP=capacities, NMASK=ST_mask, &
                      NODES=n_vertices, NARCS=n_edges, &
                      ISOURC=source_vertex, ISINK=sink_vertex, MAXFLO=i_temp, &
                      NNCUT=n_S_nodes, NACUT=n_cut_arcs, IRETN=status) @;
      IF(status!=0) THEN @;
         CALL NonCriticalError(message="Min Cut not found correctly for a subgraph", &
                               caller="SubMinCut") @;
      END IF @; 
      n_T_nodes=n_vertices-n_S_nodes @;                       
      @#endif
      
      @#if 0
         WRITE(*,*) "---" @;
         WRITE(*,*) ST_mask @;         
         WRITE(*,*) n_S_nodes, n_cut_arcs, i_temp @;
         PAUSE @;
      @#endif      

      DO edge=1, n_edges @;
         IF(heads[edge]==sink_vertex) THEN @;
            IF(ST_mask(tails[edge])) THEN @; // This cut is not valid--it contains ST arcs (of inifinite capacity)
               cut_capacity=0.0_r_wp @;
               @%% WRITE(*,*) "Invalid cut for parent_cut=", parent_cut @;
               RETURN @; // Nothing happens
            END IF @;
         ELSE IF(tails[edge]==source_vertex) THEN @;  
            IF(.NOT.(ST_mask(tails[edge]))) THEN @; // This cut is not valid--it contains ST arcs (of inifinite capacity)
               cut_capacity=0.0_r_wp @;
               @%% WRITE(*,*) "Invalid cut for parent_cut=", parent_cut @;               
               RETURN @; // Nothing happens
            END IF @;
         END IF @;
      END DO @;
                     
      cut_capacity=@E REAL(i_temp,r_wp)/normalization @; // The max-flow/min-cut value   
            
      @#if 1
      DO node=1, n_nodes @;
         vertex=nodes_renumbering[node] @;
         IF(vertex!=0) THEN @;            
            IF(.NOT.ST_mask[vertex]) THEN @; // Increment the label of the sink nodes by 1
               nodes_cuts[node]=parent_cut+1 @;
            END IF @;
         END IF @;
      END DO @;
      @#endif     
            
      _DeallocateArray(nodes_renumbering, 1_i_wp, ALLOCATED, "LabelingMinCut") @; 
      _DeallocateArray(heads, 1_i_wp, ALLOCATED, "LabelingMinCut") @;
      _DeallocateArray(tails, 1_i_wp, ALLOCATED, "LabelingMinCut") @;   
      _DeallocateArray(capacities, 1_i_wp, ALLOCATED, "LabelingMinCut") @;  
      _DeallocateArray(ST_mask, .TRUE._l_short, ALLOCATED, "LabelingMinCut") @;  
                                
   END SUBROUTINE SubMinCut @;

END SUBROUTINE LabelingMinCut @;

@%%