MODULE Graph_Algorithms
      USE Precision
      USE Error_Handling
      USE System_Monitors
      IMPLICIT NONE
      PUBLIC::ConnectedComponents
      PRIVATE
      
      CONTAINS
      
      
      SUBROUTINE ConnectedComponents(heads_tails,node_offset,labels,use_labels,&
      largest_component_label,largest_component_size)
      INTEGER(KIND=i_wp),DIMENSION(:,:),INTENT(IN)::heads_tails
      INTEGER(KIND=i_wp),INTENT(IN)::node_offset
      INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(INOUT)::labels
      LOGICAL,INTENT(IN),OPTIONAL::use_labels
      INTEGER(KIND=i_wp),INTENT(OUT),OPTIONAL::largest_component_label,largest_component_size
      
      INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::heights
      INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_arcs
      INTEGER(KIND=i_wp)::node,arc,head,tail,head_root,tail_root,&
      root,diff_heights,child,parent,&
      largest_component(1),climber,climber_parent
      INTEGER::alloc_status
      LOGICAL::reuse_labels
      
      n_special_nodes=-INT(LBOUND(labels,DIM=1),KIND=i_wp)
      n_nodes=INT(UBOUND(labels,DIM=1),KIND=i_wp)
      n_arcs=INT(SIZE(heads_tails,DIM=2),KIND=i_wp)
      
      reuse_labels=.FALSE.
      IF(PRESENT(use_labels))reuse_labels=use_labels
      IF(.NOT.reuse_labels)THEN
      DO node=-n_special_nodes,n_nodes
      labels(node)=node
      END DO
      END IF
      
      ALLOCATE(heights(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=n_nodes,mold=1_i_wp,&
      caller="ConnectedComponents",alloc_status=alloc_status)
      heights=0
      
      ForAllArcs:DO arc=1,n_arcs
      head=heads_tails(1,arc)
      tail=heads_tails(2,arc)
      
      
      climber=head
      DO WHILE(labels(climber)/=climber)
      climber=labels(climber)
      END DO
      head_root=climber
      
      climber=tail
      DO WHILE(labels(climber)/=climber)
      climber=labels(climber)
      END DO
      tail_root=climber
      
      climber=head
      DO WHILE(labels(climber)/=head_root)
      climber_parent=labels(climber)
      labels(climber)=head_root
      climber=climber_parent
      END DO
      
      climber=tail
      DO WHILE(labels(climber)/=tail_root)
      climber_parent=labels(climber)
      labels(climber)=tail_root
      climber=climber_parent
      END DO
      
      TreeUnion:IF(head_root/=tail_root)THEN
      diff_heights=(heights(head_root)-heights(tail_root))
      IF(diff_heights>0)THEN
      parent=head_root
      child=tail_root
      ELSE
      parent=tail_root
      child=head_root
      END IF
      labels(child)=parent
      IF(diff_heights==0)heights(parent)=heights(parent)+1
      END IF TreeUnion
      
      END DO ForAllArcs
      
      ForAllNodes:DO node=-n_special_nodes,n_nodes
      
      climber=node
      DO WHILE(labels(climber)/=climber)
      climber=labels(climber)
      END DO
      root=climber
      
      climber=node
      DO WHILE(labels(climber)/=root)
      climber_parent=labels(climber)
      labels(climber)=root
      climber=climber_parent
      END DO
      END DO ForAllNodes
      
      
      heights=0
      DO node=-n_special_nodes,n_nodes
      heights(labels(node))=heights(labels(node))+1
      END DO
      
      IF(PRESENT(largest_component_label).OR.PRESENT(largest_component_size))THEN
      largest_component=INT(MAXLOC(heights),KIND=i_wp)
      IF(PRESENT(largest_component_label))&
      largest_component_label=labels(largest_component(1)-1-node_offset)
      IF(PRESENT(largest_component_size))&
      largest_component_size=heights(largest_component(1)-1-node_offset)
      END IF
      
      DEALLOCATE(heights)
      CALL RecordAllocation(n_elements=-n_nodes,mold=1_i_wp,caller="ConnectedComponents")
      
      END SUBROUTINE ConnectedComponents
      
      
      
      END MODULE Graph_Algorithms