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
      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
      
      
      SUBROUTINE DijkstraShortestPaths(node_offset,heads_tails,special_heads_tails,arcs_lengths,&
      source_sink,shortest_paths,n_paths,shortest_lengths,queue_type)
      
      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
      INTEGER(KIND=i_wp),DIMENSION(2),INTENT(IN)::source_sink
      INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(OUT)::shortest_paths
      INTEGER(KIND=i_byte),INTENT(IN)::n_paths
      REAL(KIND=r_wp),DIMENSION(:),INTENT(OUT),OPTIONAL::shortest_lengths
      CHARACTER,INTENT(IN),OPTIONAL::queue_type
      
      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
      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'
      IF(PRESENT(queue_type))queue=queue_type
      
      n_special_nodes=-INT(LBOUND(shortest_paths,DIM=1),KIND=i_wp)
      n_nodes=INT(UBOUND(shortest_paths,DIM=1),KIND=i_wp)
      n_special_arcs=INT(SIZE(special_heads_tails,DIM=2),KIND=i_wp)
      n_arcs=INT(SIZE(heads_tails,DIM=2),KIND=i_wp)
      
      IF(ANY(source_sink>=0))THEN
      
      CALL NonCriticalError("Source or sink is not a special node","DijkstraShortestPaths")
      RETURN
      END IF
      n_vertices=n_nodes+2
      root=n_nodes+1
      n_edges=2*(n_arcs+n_special_arcs)
      
      
      IF(.NOT.ALLOCATED(my_neighbours))THEN
      ALLOCATE(my_neighbours(0:n_vertices+1),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_vertices+1)-(0)+1,mold=1_i_wp,&
      caller="DijkstraShortestPaths",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ALLOCATED(distances))THEN
      ALLOCATE(distances(0:n_vertices+1),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_vertices+1)-(0)+1,mold=1_i_wp,&
      caller="DijkstraShortestPaths",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ALLOCATED(parents))THEN
      ALLOCATE(parents(0:n_vertices+1),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_vertices+1)-(0)+1,mold=1_i_wp,&
      caller="DijkstraShortestPaths",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ALLOCATED(neighbours))THEN
      ALLOCATE(neighbours(0:n_edges),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_edges)-(0)+1,mold=1_i_wp,&
      caller="DijkstraShortestPaths",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ALLOCATED(incident_arcs))THEN
      ALLOCATE(incident_arcs(0:n_edges),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_edges)-(0)+1,mold=1_i_wp,&
      caller="DijkstraShortestPaths",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ALLOCATED(edge_lengths))THEN
      ALLOCATE(edge_lengths(0:n_edges),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_edges)-(0)+1,mold=1_i_wp,&
      caller="DijkstraShortestPaths",alloc_status=alloc_status)
      END IF
      
      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)
      max_edge_weight=HUGE(1_i_wp)/(10*n_vertices)
      DO edge=1,n_edges
      arc=incident_arcs(edge)
      edge_lengths(edge)=INT(REAL(max_edge_weight)*arcs_lengths(arc)/max_weight,i_wp)+1_i_wp
      END DO
      
      ST_degrees=0
      STArcs:DO s_or_t=1,2
      my_neighbours(n_nodes+s_or_t)=n_edges+1
      DO special_arc=1,n_special_arcs
      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
      
      
      
      IF(.NOT.ALLOCATED(temp_array_1))THEN
      ALLOCATE(temp_array_1(0:n_vertices+1),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_vertices+1)-(0)+1,mold=1_i_wp,&
      caller="DijkstraShortestPaths",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ALLOCATED(temp_array_2))THEN
      ALLOCATE(temp_array_2(0:n_vertices+1),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_vertices+1)-(0)+1,mold=1_i_wp,&
      caller="DijkstraShortestPaths",alloc_status=alloc_status)
      END IF
      
      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)
      
      
      IF(ALLOCATED(temp_array_1))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(temp_array_1),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(temp_array_1,STAT=alloc_status)
      END IF
      
      IF(ALLOCATED(temp_array_2))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(temp_array_2),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(temp_array_2,STAT=alloc_status)
      END IF
      
      
      IF(.NOT.ALLOCATED(paths_lengths))THEN
      ALLOCATE(paths_lengths(1:ST_degrees(2)),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(ST_degrees(2))-(1)+1,mold=1_i_wp,&
      caller="DijkstraShortestPaths",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ALLOCATED(paths_ordering))THEN
      ALLOCATE(paths_ordering(1:ST_degrees(2)),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(ST_degrees(2))-(1)+1,mold=1_i_wp,&
      caller="DijkstraShortestPaths",alloc_status=alloc_status)
      END IF
      
      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))
      
      
      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
      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)=REAL(min_distance)/REAL(distances(vertex))
      MarkPath:DO
      counter=counter+1
      node=vertex
      IF((shortest_paths(node)/=0).OR.(vertex==source_vertex))EXIT MarkPath
      shortest_paths(node)=path
      vertex=parents(vertex)
      END DO MarkPath
      END DO SelectedPaths
      
      IF(ALLOCATED(paths_lengths))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(paths_lengths),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(paths_lengths,STAT=alloc_status)
      END IF
      
      IF(ALLOCATED(paths_ordering))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(paths_ordering),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(paths_ordering,STAT=alloc_status)
      END IF
      
      
      IF(ALLOCATED(neighbours))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(neighbours),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(neighbours,STAT=alloc_status)
      END IF
      
      IF(ALLOCATED(parents))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(parents),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(parents,STAT=alloc_status)
      END IF
      
      IF(ALLOCATED(distances))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(distances),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(distances,STAT=alloc_status)
      END IF
      
      
      IF(ALLOCATED(my_neighbours))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(my_neighbours),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(my_neighbours,STAT=alloc_status)
      END IF
      
      IF(ALLOCATED(incident_arcs))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(incident_arcs),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(incident_arcs,STAT=alloc_status)
      END IF
      
      IF(ALLOCATED(edge_lengths))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(edge_lengths),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(edge_lengths,STAT=alloc_status)
      END IF
      
      END SUBROUTINE DijkstraShortestPaths
      
      
      
      SUBROUTINE LabelingMinCut(node_offset,arc_offset,heads_tails,arcs_capacities,&
      source_sink,nodes_cuts,arcs_cuts,n_cuts,cuts_capacities)
      
      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
      INTEGER(KIND=i_wp),DIMENSION(2),INTENT(IN)::source_sink
      INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(OUT)::nodes_cuts
      INTEGER(KIND=i_byte),DIMENSION(:),INTENT(OUT)::arcs_cuts
      INTEGER(KIND=i_byte),INTENT(IN)::n_cuts
      REAL(KIND=r_wp),DIMENSION(:),INTENT(OUT),OPTIONAL::cuts_capacities
      
      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
      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
      n_special_nodes=-INT(LBOUND(nodes_cuts,DIM=1),KIND=i_wp)
      n_nodes=INT(UBOUND(nodes_cuts,DIM=1),KIND=i_wp)
      n_special_arcs=-INT(LBOUND(heads_tails,DIM=2),KIND=i_wp)
      n_arcs=INT(UBOUND(heads_tails,DIM=2),KIND=i_wp)
      IF(ANY(source_sink>=0))THEN
      
      CALL NonCriticalError("Source or sink is not a special node","LabelingMinCut")
      RETURN
      END IF
      
      max_capacity=MAXVAL(arcs_capacities)
      normalization=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
      nodes_cuts(source_sink(1))=-1
      nodes_cuts(source_sink(2))=HUGE(1_i_byte)
      arcs_cuts=0
      
      cut_level=0
      low_cut=0
      high_cut=0
      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)
      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
      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
      END IF
      END DO MarkArcs
      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
      
      
      
      IF(.NOT.ALLOCATED(cuts_reordering))THEN
      ALLOCATE(cuts_reordering(1:n_cuts),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_cuts)-(1)+1,mold=1_i_wp,&
      caller="LabelingMinCut",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ALLOCATED(cuts_renumbering))THEN
      ALLOCATE(cuts_renumbering(1:n_cuts),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_cuts)-(1)+1,mold=1_i_byte,&
      caller="LabelingMinCut",alloc_status=alloc_status)
      END IF
      CALL QuickRank(array=cuts_capacities(1:n_cuts),permutation=cuts_reordering)
      
      DO cut=1,n_cuts
      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))
      
      
      DO arc=1,n_arcs
      IF(arcs_cuts(arc)>0)THEN
      arcs_cuts(arc)=cuts_renumbering(arcs_cuts(arc))
      END IF
      END DO
      
      
      IF(ALLOCATED(cuts_reordering))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(cuts_reordering),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(cuts_reordering,STAT=alloc_status)
      END IF
      
      IF(ALLOCATED(cuts_renumbering))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(cuts_renumbering),KIND=i_wp),mold=1_i_byte)
      DEALLOCATE(cuts_renumbering,STAT=alloc_status)
      END IF
      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
      REAL(KIND=r_wp),INTENT(IN),OPTIONAL::capacity_normalization
      INTEGER(KIND=i_byte),INTENT(IN)::parent_cut
      INTEGER(KIND=i_byte),DIMENSION(-node_offset:),INTENT(INOUT)::cuts_labels
      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
      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=-INT(LBOUND(cuts_labels,DIM=1),KIND=i_wp)
      n_nodes=INT(UBOUND(cuts_labels,DIM=1),KIND=i_wp)
      n_special_arcs=-INT(LBOUND(heads_tails,DIM=2),KIND=i_wp)
      n_arcs=INT(UBOUND(heads_tails,DIM=2),KIND=i_wp)
      
      
      IF(.NOT.ALLOCATED(nodes_renumbering))THEN
      ALLOCATE(nodes_renumbering(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_wp,&
      caller="LabelingMinCut",alloc_status=alloc_status)
      END IF
      
      nodes_renumbering=0
      counter=0
      DO node=-n_special_nodes,n_nodes
      IF(cuts_labels(node)==parent_cut)THEN
      counter=counter+1
      nodes_renumbering(node)=counter
      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
      head=heads_tails(1,arc)
      tail=heads_tails(2,arc)
      
      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
      
      
      IF(.NOT.ALLOCATED(heads))THEN
      ALLOCATE(heads(1:n_edges),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_edges)-(1)+1,mold=1_i_wp,&
      caller="LabelingMinCut",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ALLOCATED(tails))THEN
      ALLOCATE(tails(1:n_edges),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_edges)-(1)+1,mold=1_i_wp,&
      caller="LabelingMinCut",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ALLOCATED(capacities))THEN
      ALLOCATE(capacities(1:n_edges),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_edges)-(1)+1,mold=1_i_wp,&
      caller="LabelingMinCut",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ALLOCATED(ST_mask))THEN
      ALLOCATE(ST_mask(1:n_vertices),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_vertices)-(1)+1,mold=.TRUE._l_short,&
      caller="LabelingMinCut",alloc_status=alloc_status)
      END IF
      
      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
      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
      END IF
      
      ELSE IF(cuts_labels(head)==parent_cut)THEN
      
      edge=edge+1
      IF(cuts_labels(tail)<parent_cut)THEN
      n_S_arcs=n_S_arcs+1
      tails(edge)=source_vertex
      heads(edge)=nodes_renumbering(head)
      ELSE
      n_T_arcs=n_T_arcs+1
      heads(edge)=sink_vertex
      tails(edge)=nodes_renumbering(head)
      END IF
      capacities(edge)=inf_capacity
      
      ELSE IF(cuts_labels(tail)==parent_cut)THEN
      
      edge=edge+1
      IF(cuts_labels(head)<parent_cut)THEN
      n_S_arcs=n_S_arcs+1
      tails(edge)=source_vertex
      heads(edge)=nodes_renumbering(tail)
      ELSE
      n_T_arcs=n_T_arcs+1
      heads(edge)=sink_vertex
      tails(edge)=nodes_renumbering(tail)
      END IF
      capacities(edge)=inf_capacity
      
      END IF
      
      END DO
      
      
      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
      
      
      DO edge=1,n_edges
      IF(heads(edge)==sink_vertex)THEN
      IF(ST_mask(tails(edge)))THEN
      cut_capacity=0.0_r_wp
      RETURN
      END IF
      ELSE IF(tails(edge)==source_vertex)THEN
      IF(.NOT.(ST_mask(tails(edge))))THEN
      cut_capacity=0.0_r_wp
      RETURN
      END IF
      END IF
      END DO
      
      cut_capacity=REAL(i_temp,r_wp)/normalization
      
      DO node=1,n_nodes
      vertex=nodes_renumbering(node)
      IF(vertex/=0)THEN
      IF(.NOT.ST_mask(vertex))THEN
      nodes_cuts(node)=parent_cut+1
      END IF
      END IF
      END DO
      
      
      IF(ALLOCATED(nodes_renumbering))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_renumbering),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(nodes_renumbering,STAT=alloc_status)
      END IF
      
      IF(ALLOCATED(heads))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(heads),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(heads,STAT=alloc_status)
      END IF
      
      IF(ALLOCATED(tails))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(tails),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(tails,STAT=alloc_status)
      END IF
      
      IF(ALLOCATED(capacities))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(capacities),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(capacities,STAT=alloc_status)
      END IF
      
      IF(ALLOCATED(ST_mask))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(ST_mask),KIND=i_wp),mold=.TRUE._l_short)
      DEALLOCATE(ST_mask,STAT=alloc_status)
      END IF
      
      END SUBROUTINE SubMinCut
      
      END SUBROUTINE LabelingMinCut
      
      
      
      END MODULE Graph_Algorithms_F77