MODULE SCOTCH_Constants
      USE Precision
      IMPLICIT NONE
      PUBLIC
      
      INTEGER,PARAMETER::i_scotch=i_word,SCOTCH_MAX_LEN=200
      
      INTEGER,PARAMETER::SCOTCH_ARCHDIM=7,SCOTCH_GEOMDIM=15,SCOTCH_GRAPHDIM=14,&
      SCOTCH_MAPDIM=32,SCOTCH_ORDERDIM=9,SCOTCH_STRATDIM=1
      
      
      END MODULE SCOTCH_Constants
      
      MODULE SCOTCH_Interface
      USE Precision
      USE SCOTCH_Constants
      USE Error_Handling
      USE System_Monitors
      USE Graph_Algorithms
      USE Network_Data_Types
      IMPLICIT NONE
      PRIVATE
      PUBLIC::SCOTCH_InitializeMappingOrdering,SCOTCH_ComputeMappingOrdering,&
      SCOTCH_DestroyMappingOrdering
      
      ML_EXTERNAL::SCOTCH_graphInit,SCOTCH_graphExit,SCOTCH_graphWrite,&
      SCOTCH_graphBuild,SCOTCH_graphMap,SCOTCH_graphOrder,&
      SCOTCH_stratInit,SCOTCH_stratExit,SCOTCH_stratMap,SCOTCH_stratOrder,&
      SCOTCH_archInit,SCOTCH_archExit,SCOTCH_archAlgorithmic
      INTEGER::SCOTCH_graphInit,SCOTCH_graphBuild,SCOTCH_graphMap,SCOTCH_graphOrder,&
      SCOTCH_stratInit,SCOTCH_stratMap,SCOTCH_stratOrder,&
      SCOTCH_archInit,SCOTCH_archAlgorithmic
      
      PUBLIC::SCOTCH_Graph_,SCOTCH_Mapping_,&
      SCOTCH_Ordering_,SCOTCH_Mapping_Ordering
      
      TYPE SCOTCH_Graph_
      PRIVATE
      INTEGER(KIND=i_scotch),DIMENSION(SCOTCH_GRAPHDIM)::graph
      INTEGER(KIND=i_scotch)::n_vertices=0,n_edges=0
      INTEGER(KIND=i_scotch),DIMENSION(:),POINTER::verttab,edgetab,velotab,&
      edlotab,incitab
      ENDTYPE
      TYPE SCOTCH_Mapping_
      INTEGER(KIND=i_scotch),DIMENSION(SCOTCH_MAPDIM)::mapping
      INTEGER(KIND=i_scotch),DIMENSION(SCOTCH_STRATDIM)::strategy
      INTEGER(KIND=i_scotch),DIMENSION(SCOTCH_ARCHDIM)::architecture
      CHARACTER(LEN=SCOTCH_MAX_LEN)::architecture_string=""
      
      CHARACTER(LEN=SCOTCH_MAX_LEN)::mapping_string=""
      ENDTYPE
      TYPE SCOTCH_Ordering_
      INTEGER(KIND=i_scotch),DIMENSION(SCOTCH_ORDERDIM)::ordering
      INTEGER(KIND=i_scotch),DIMENSION(SCOTCH_STRATDIM)::strategy
      INTEGER(KIND=i_scotch)::n_blocks=0
      INTEGER(KIND=i_scotch),DIMENSION(:),POINTER::block_pointers
      
      CHARACTER(LEN=SCOTCH_MAX_LEN)::ordering_string=""
      ENDTYPE
      
      TYPE SCOTCH_Mapping_Ordering
      LOGICAL::map_graph=.TRUE.,order_graph=.FALSE.
      TYPE(Directed_Graph),POINTER::graph
      TYPE(Support_Tree_Mapping),POINTER::ST_mapping
      TYPE(Fill_Reducing_Ordering),POINTER::graph_ordering
      TYPE(SCOTCH_Graph_)::scotch_graph
      TYPE(SCOTCH_Mapping_)::scotch_mapping
      TYPE(SCOTCH_Ordering_)::scotch_ordering
      LOGICAL::allocated_mapping=.FALSE.,allocated_ordering=.FALSE.
      ENDTYPE
      
      CONTAINS
      
      
      SUBROUTINE SCOTCH_InitializeMappingOrdering(mapping_ordering)
      IMPLICIT NONE
      TYPE(SCOTCH_Mapping_Ordering),INTENT(INOUT),TARGET::mapping_ordering
      
      TYPE(Directed_Graph),POINTER::graph
      INTEGER::alloc_status,result
      INTEGER(KIND=i_scotch)::n_vertices,n_edges
      INTEGER(KIND=i_wp)::n_special_nodes,n_nodes,n_special_arcs,n_arcs
      LOGICAL::valid_graph
      
      graph=>mapping_ordering%graph
      
      valid_graph=.TRUE.
      IF(.NOT.ASSOCIATED(graph))THEN
      valid_graph=.FALSE.
      ELSE
      IF(.NOT.ASSOCIATED(graph%heads_tails))valid_graph=.FALSE.
      END IF
      IF(.NOT.valid_graph)THEN
      CALL CriticalError(message="Incomplete or no graph passed to SCOTCH library",&
      caller="SCOTCH_InitializeMappingOrdering")
      RETURN
      END IF
      
      n_special_nodes=graph%n_special_nodes
      n_nodes=graph%n_nodes
      n_special_arcs=graph%n_special_arcs
      n_arcs=graph%n_arcs
      
      n_vertices=n_nodes+n_special_nodes+1
      n_edges=2*(n_arcs+n_special_arcs+1)
      
      result=SCOTCH_graphInit(mapping_ordering%scotch_graph%graph)
      result=SCOTCH_archInit(mapping_ordering%scotch_mapping%architecture)
      IF(mapping_ordering%map_graph)result=SCOTCH_stratInit(mapping_ordering%scotch_mapping%strategy)
      IF(mapping_ordering%order_graph)result=SCOTCH_stratInit(mapping_ordering%scotch_ordering%strategy)
      
      
      IF(.NOT.ASSOCIATED(mapping_ordering%scotch_graph%verttab))THEN
      ALLOCATE(mapping_ordering%scotch_graph%verttab(0:n_vertices+1),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_vertices+1)-(0)+1,mold=1_i_scotch,&
      caller="SCOTCH_InitializeMappingOrdering",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ASSOCIATED(mapping_ordering%scotch_graph%edgetab))THEN
      ALLOCATE(mapping_ordering%scotch_graph%edgetab(0:n_edges),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_edges)-(0)+1,mold=1_i_scotch,&
      caller="SCOTCH_InitializeMappingOrdering",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ASSOCIATED(mapping_ordering%scotch_graph%incitab))THEN
      ALLOCATE(mapping_ordering%scotch_graph%incitab(0:n_edges),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_edges)-(0)+1,mold=1_i_scotch,&
      caller="SCOTCH_InitializeMappingOrdering",alloc_status=alloc_status)
      END IF
      
      CALL CreateAdjacencyArrays(node_offset=n_special_nodes,arc_offset=n_special_arcs,&
      heads_tails=graph%heads_tails,neighbours=mapping_ordering%scotch_graph%edgetab,&
      my_neighbours=mapping_ordering%scotch_graph%verttab,incident_arcs=mapping_ordering%scotch_graph%incitab,&
      n_edges=n_edges,include_loops=.FALSE.)
      mapping_ordering%scotch_graph%n_vertices=n_vertices
      mapping_ordering%scotch_graph%n_edges=n_edges
      
      IF(mapping_ordering%map_graph)THEN
      IF(.NOT.ASSOCIATED(mapping_ordering%ST_mapping))THEN
      ALLOCATE(mapping_ordering%ST_mapping,STAT=alloc_status)
      mapping_ordering%allocated_mapping=.TRUE.
      END IF
      
      IF(.NOT.ASSOCIATED(mapping_ordering%scotch_graph%edlotab))THEN
      ALLOCATE(mapping_ordering%scotch_graph%edlotab(0:n_edges),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_edges)-(0)+1,mold=1_i_scotch,&
      caller="SCOTCH_InitializeMappingOrdering",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(mapping_ordering%ST_mapping%arcs_weights))THEN
      ALLOCATE(mapping_ordering%ST_mapping%arcs_weights(-n_special_arcs:n_arcs),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1.0_r_wp,&
      caller="SCOTCH_InitializeMappingOrdering",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(mapping_ordering%ST_mapping%nodes_mapping))THEN
      ALLOCATE(mapping_ordering%ST_mapping%nodes_mapping(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_wp,&
      caller="SCOTCH_InitializeMappingOrdering",alloc_status=alloc_status)
      END IF
      END IF
      
      IF(mapping_ordering%order_graph)THEN
      IF(.NOT.ASSOCIATED(mapping_ordering%graph_ordering))THEN
      ALLOCATE(mapping_ordering%graph_ordering,STAT=alloc_status)
      mapping_ordering%allocated_ordering=.TRUE.
      END IF
      
      IF(.NOT.ASSOCIATED(mapping_ordering%scotch_ordering%block_pointers))THEN
      ALLOCATE(mapping_ordering%scotch_ordering%block_pointers(0:n_vertices+1),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_vertices+1)-(0)+1,mold=1_i_scotch,&
      caller="SCOTCH_InitializeMappingOrdering",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(mapping_ordering%graph_ordering%nodes_reordering))THEN
      ALLOCATE(mapping_ordering%graph_ordering%nodes_reordering(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_wp,&
      caller="SCOTCH_InitializeMappingOrdering",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(mapping_ordering%graph_ordering%nodes_renumbering))THEN
      ALLOCATE(mapping_ordering%graph_ordering%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="SCOTCH_InitializeMappingOrdering",alloc_status=alloc_status)
      END IF
      END IF
      
      END SUBROUTINE SCOTCH_InitializeMappingOrdering
      
      
      
      SUBROUTINE SCOTCH_ComputeMappingOrdering(mapping_ordering,mapping_offset)
      IMPLICIT NONE
      TYPE(SCOTCH_Mapping_Ordering),INTENT(INOUT),TARGET::mapping_ordering
      INTEGER(KIND=i_wp),INTENT(IN),OPTIONAL::mapping_offset
      
      TYPE(Directed_Graph),POINTER::graph
      INTEGER(KIND=i_wp)::arc
      INTEGER(KIND=i_scotch)::edge,max_scotch_weight
      REAL(KIND=r_wp)::max_weight
      INTEGER::alloc_status,result
      
      graph=>mapping_ordering%graph
      
      MapGraph:IF(mapping_ordering%map_graph)THEN
      
      result=SCOTCH_archAlgorithmic(mapping_ordering%scotch_mapping%architecture,mapping_ordering%scotch_mapping%architecture_stri&
     &ng//CHAR(0))
      result=SCOTCH_stratMap(mapping_ordering%scotch_mapping%strategy,mapping_ordering%scotch_mapping%mapping_string//CHAR(0))
      max_weight=MAXVAL(mapping_ordering%ST_mapping%arcs_weights)
      max_scotch_weight=HUGE(1_i_scotch)/10/(mapping_ordering%scotch_graph%n_vertices+mapping_ordering%scotch_graph%n_edges)
      DO edge=1,mapping_ordering%scotch_graph%n_edges
      arc=mapping_ordering%scotch_graph%incitab(edge)
      mapping_ordering%scotch_graph%edlotab(edge)=INT(REAL(max_scotch_weight)*mapping_ordering%ST_mapping%arcs_weights(arc)/max_we&
     &ight,&
      i_scotch)+1_i_scotch
      END DO
      
      
      result=SCOTCH_graphBuild(mapping_ordering%scotch_graph%graph,%VAL(INT(mapping_ordering%scotch_graph%n_vertices,i_scotch)),&
      mapping_ordering%scotch_graph%verttab(1),%VAL(0),%VAL(0),%VAL(INT(mapping_ordering%scotch_graph%n_edges,i_scotch)),&
      mapping_ordering%scotch_graph%edgetab(1),mapping_ordering%scotch_graph%edlotab(1),%VAL(1_i_scotch),1_i_scotch)
      CALL SCOTCH_graphWrite(mapping_ordering%scotch_graph%graph,"graph.out"//CHAR(0))
      result=SCOTCH_graphMap(mapping_ordering%scotch_graph%graph,mapping_ordering%scotch_mapping%strategy,&
      mapping_ordering%scotch_mapping%architecture,mapping_ordering%ST_mapping%nodes_mapping(-graph%n_special_nodes))
      IF(PRESENT(mapping_offset))THEN
      mapping_ordering%ST_mapping%nodes_mapping=mapping_ordering%ST_mapping%nodes_mapping+mapping_offset
      END IF
      
      END IF MapGraph
      
      OrderGraph:IF(mapping_ordering%order_graph)THEN
      
      result=SCOTCH_stratOrder(mapping_ordering%scotch_ordering%strategy,mapping_ordering%scotch_ordering%ordering_string//CHAR(0)&
     &)
      result=SCOTCH_graphBuild(mapping_ordering%scotch_graph%graph,%VAL(INT(mapping_ordering%scotch_graph%n_vertices,i_scotch)),&
      mapping_ordering%scotch_graph%verttab(1),%VAL(0),%VAL(0),%VAL(INT(mapping_ordering%scotch_graph%n_edges,i_scotch)),&
      mapping_ordering%scotch_graph%edgetab(1),%VAL(0),%VAL(1_i_scotch),3_i_scotch)
      result=SCOTCH_graphOrder(mapping_ordering%scotch_graph%graph,mapping_ordering%scotch_ordering%strategy,&
      mapping_ordering%scotch_ordering%n_blocks,mapping_ordering%scotch_ordering%block_pointers(1),&
      mapping_ordering%graph_ordering%nodes_renumbering(-graph%n_special_nodes),mapping_ordering%graph_ordering%nodes_reordering(-&
     &graph%n_special_nodes))
      mapping_ordering%graph_ordering%nodes_reordering=mapping_ordering%graph_ordering%nodes_reordering-graph%n_special_nodes-1
      mapping_ordering%graph_ordering%nodes_renumbering=mapping_ordering%graph_ordering%nodes_renumbering-graph%n_special_nodes-1
      END IF OrderGraph
      
      END SUBROUTINE SCOTCH_ComputeMappingOrdering
      
      
      
      SUBROUTINE SCOTCH_DestroyMappingOrdering(mapping_ordering)
      IMPLICIT NONE
      TYPE(SCOTCH_Mapping_Ordering),INTENT(INOUT),TARGET::mapping_ordering
      
      INTEGER::alloc_status,result
      
      CALL SCOTCH_graphExit(mapping_ordering%scotch_graph%graph)
      CALL SCOTCH_archExit(mapping_ordering%scotch_mapping%architecture)
      IF(mapping_ordering%map_graph)CALL SCOTCH_stratExit(mapping_ordering%scotch_mapping%strategy)
      IF(mapping_ordering%order_graph)CALL SCOTCH_stratExit(mapping_ordering%scotch_ordering%strategy)
      
      IF(mapping_ordering%allocated_mapping)THEN
      
      IF(ASSOCIATED(mapping_ordering%ST_mapping%arcs_weights))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(mapping_ordering%ST_mapping%arcs_weights),KIND=i_wp),mold=1.0_r_wp)
      DEALLOCATE(mapping_ordering%ST_mapping%arcs_weights,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(mapping_ordering%ST_mapping%nodes_mapping))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(mapping_ordering%ST_mapping%nodes_mapping),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(mapping_ordering%ST_mapping%nodes_mapping,STAT=alloc_status)
      END IF
      DEALLOCATE(mapping_ordering%ST_mapping,STAT=alloc_status)
      END IF
      
      IF(mapping_ordering%allocated_ordering)THEN
      
      IF(ASSOCIATED(mapping_ordering%graph_ordering%nodes_reordering))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(mapping_ordering%graph_ordering%nodes_reordering),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(mapping_ordering%graph_ordering%nodes_reordering,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(mapping_ordering%graph_ordering%nodes_renumbering))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(mapping_ordering%graph_ordering%nodes_renumbering),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(mapping_ordering%graph_ordering%nodes_renumbering,STAT=alloc_status)
      END IF
      DEALLOCATE(mapping_ordering%graph_ordering,STAT=alloc_status)
      END IF
      
      
      IF(ASSOCIATED(mapping_ordering%scotch_graph%verttab))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(mapping_ordering%scotch_graph%verttab),KIND=i_wp),mold=1_i_scotch)
      DEALLOCATE(mapping_ordering%scotch_graph%verttab,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(mapping_ordering%scotch_graph%edgetab))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(mapping_ordering%scotch_graph%edgetab),KIND=i_wp),mold=1_i_scotch)
      DEALLOCATE(mapping_ordering%scotch_graph%edgetab,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(mapping_ordering%scotch_graph%incitab))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(mapping_ordering%scotch_graph%incitab),KIND=i_wp),mold=1_i_scotch)
      DEALLOCATE(mapping_ordering%scotch_graph%incitab,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(mapping_ordering%scotch_graph%edlotab))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(mapping_ordering%scotch_graph%edlotab),KIND=i_wp),mold=1_i_scotch)
      DEALLOCATE(mapping_ordering%scotch_graph%edlotab,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(mapping_ordering%scotch_ordering%block_pointers))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(mapping_ordering%scotch_ordering%block_pointers),KIND=i_wp),mold=1_i_scotch)
      DEALLOCATE(mapping_ordering%scotch_ordering%block_pointers,STAT=alloc_status)
      END IF
      
      END SUBROUTINE SCOTCH_DestroyMappingOrdering
      
      
      
      END MODULE SCOTCH_Interface