LCOV - code coverage report
Current view: top level - src - kg_vertex_coloring_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 82.5 % 251 207
Test Date: 2025-07-25 12:55:17 Functions: 73.3 % 15 11

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Routines for a Kim-Gordon-like partitioning into molecular subunits
      10              : !>        unsing a vertex coloring algorithm (DSATUR) to find non-interating
      11              : !>        subsets, such that two molecules within the same subset have
      12              : !>        small/zero overlap (in other words: this molecular pair is not present in
      13              : !>        the neighborlist sab_orb for the current value of EPS_DEFAULT)
      14              : !> \par History
      15              : !>       2012.07 created [Martin Haeufel]
      16              : !>         2013.11 Added pair switching and revised Dsatur [Samuel Andermatt]
      17              : !> \author Martin Haeufel
      18              : ! **************************************************************************************************
      19              : MODULE kg_vertex_coloring_methods
      20              :    USE bibliography,                    ONLY: Andermatt2016,&
      21              :                                               Brelaz1979,&
      22              :                                               cite_reference
      23              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      24              :                                               cp_logger_get_default_unit_nr,&
      25              :                                               cp_logger_type
      26              :    USE cp_min_heap,                     ONLY: cp_heap_fill,&
      27              :                                               cp_heap_new,&
      28              :                                               cp_heap_pop,&
      29              :                                               cp_heap_release,&
      30              :                                               cp_heap_reset_node,&
      31              :                                               cp_heap_type,&
      32              :                                               valt
      33              :    USE input_constants,                 ONLY: kg_color_dsatur,&
      34              :                                               kg_color_greedy
      35              :    USE kg_environment_types,            ONLY: kg_environment_type
      36              :    USE kinds,                           ONLY: int_4,&
      37              :                                               int_8
      38              : #include "./base/base_uses.f90"
      39              : 
      40              :    IMPLICIT NONE
      41              : 
      42              :    PRIVATE
      43              : 
      44              :    TYPE vertex
      45              :       INTEGER :: id = -1
      46              :       INTEGER :: color = -1
      47              :       INTEGER :: degree = -1 ! degree (number of neighbors)
      48              :       INTEGER :: dsat = -1 ! degree of saturation
      49              :       LOGICAL, ALLOCATABLE, DIMENSION(:)       :: color_present ! the colors present on neighbors
      50              :       TYPE(vertex_p_type), DIMENSION(:), POINTER :: neighbors => NULL()
      51              :    END TYPE vertex
      52              : 
      53              :    TYPE vertex_p_type
      54              :       TYPE(vertex), POINTER :: vertex => NULL()
      55              :    END TYPE vertex_p_type
      56              : 
      57              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kg_vertex_coloring_methods'
      58              : 
      59              :    PUBLIC :: kg_vertex_coloring
      60              : 
      61              : CONTAINS
      62              : ! **************************************************************************************************
      63              : !> \brief ...
      64              : !> \param kg_env ...
      65              : !> \param pairs ...
      66              : !> \param graph ...
      67              : ! **************************************************************************************************
      68           41 :    SUBROUTINE kg_create_graph(kg_env, pairs, graph)
      69              :       TYPE(kg_environment_type), POINTER                 :: kg_env
      70              :       INTEGER(KIND=int_4), ALLOCATABLE, &
      71              :          DIMENSION(:, :), INTENT(IN)                     :: pairs
      72              :       TYPE(vertex_p_type), DIMENSION(:), POINTER         :: graph
      73              : 
      74              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'kg_create_graph'
      75              : 
      76              :       INTEGER                                            :: handle, i, imol, ineighbor, jmol, kmol, &
      77              :                                                             maxdegree, nneighbors, nnodes
      78              : 
      79           41 :       CALL timeset(routineN, handle)
      80              : 
      81           41 :       CALL cite_reference(Andermatt2016)
      82              : 
      83              : ! The array pairs contains all interacting (overlapping) pairs of molecules.
      84              : ! It is assumed to be ordered in the following way:
      85              : ! (1,2), (1,3), (1,4), ..., (2,3), (2,4), ...
      86              : ! There are no entries (i,i)
      87              : ! get the number of nodes = total number of molecules
      88              : 
      89           41 :       nnodes = SIZE(kg_env%molecule_set)
      90              : 
      91           41 :       NULLIFY (graph)
      92          216 :       ALLOCATE (graph(nnodes))
      93              : 
      94              :       ! allocate and initialize all vertices
      95          134 :       DO i = 1, nnodes
      96           93 :          ALLOCATE (graph(i)%vertex)
      97           93 :          graph(i)%vertex%id = i ! id = imol (molecular index)
      98           93 :          graph(i)%vertex%color = 0 ! means vertex is not colored yet
      99           93 :          graph(i)%vertex%dsat = 0 ! no colored neighbors yet
     100          134 :          graph(i)%vertex%degree = 0 ! as needed for maxdegree....
     101              :       END DO
     102              : 
     103              :       ! allocate the neighbor lists
     104           41 :       imol = 0
     105              : 
     106           41 :       maxdegree = 0
     107              : 
     108          143 :       DO i = 1, SIZE(pairs, 2)
     109          102 :          jmol = pairs(1, i)
     110              :          ! counting loop
     111          102 :          IF (jmol .NE. imol) THEN
     112           81 :             IF (imol .NE. 0) THEN
     113          190 :                ALLOCATE (graph(imol)%vertex%neighbors(nneighbors))
     114           44 :                graph(imol)%vertex%degree = nneighbors
     115           44 :                IF (nneighbors .GT. maxdegree) maxdegree = nneighbors
     116              :             END IF
     117              :             imol = jmol
     118              :             nneighbors = 0
     119              :          END IF
     120          143 :          nneighbors = nneighbors + 1
     121              :       END DO
     122              : 
     123           41 :       IF (imol .NE. 0) THEN
     124          155 :          ALLOCATE (graph(imol)%vertex%neighbors(nneighbors))
     125           37 :          graph(imol)%vertex%degree = nneighbors
     126           37 :          IF (nneighbors .GT. maxdegree) maxdegree = nneighbors
     127              :       END IF
     128              : 
     129           41 :       kg_env%maxdegree = maxdegree
     130              : 
     131              :       ! there can be now some nodes that have no neighbors, thus vertex%neighbors
     132              :       ! is NOT ASSOCIATED
     133              : 
     134              :       ! now add the neighbors - if there are any
     135           41 :       imol = 0
     136           41 :       ineighbor = 0
     137              : 
     138          143 :       DO i = 1, SIZE(pairs, 2)
     139          102 :          jmol = pairs(1, i)
     140          102 :          IF (jmol .NE. imol) THEN
     141           81 :             ineighbor = 0
     142           81 :             imol = jmol
     143              :          END IF
     144          102 :          ineighbor = ineighbor + 1
     145          102 :          kmol = pairs(2, i)
     146          143 :          graph(imol)%vertex%neighbors(ineighbor)%vertex => graph(kmol)%vertex
     147              :       END DO
     148              : 
     149          134 :       DO i = 1, SIZE(graph)
     150          134 :          IF (graph(i)%vertex%degree > 0) THEN
     151           81 :             ALLOCATE (graph(i)%vertex%color_present(100))
     152         8181 :             graph(i)%vertex%color_present(:) = .FALSE.
     153              :          END IF
     154              :       END DO
     155              : 
     156           41 :       CALL timestop(handle)
     157              : 
     158           41 :    END SUBROUTINE
     159              : 
     160              :    ! greedy algorithm
     161              : ! **************************************************************************************************
     162              : !> \brief ...
     163              : !> \param graph ...
     164              : !> \param maxdegree ...
     165              : !> \param ncolors ...
     166              : ! **************************************************************************************************
     167            4 :    SUBROUTINE color_graph_greedy(graph, maxdegree, ncolors)
     168              :       TYPE(vertex_p_type), DIMENSION(:), POINTER         :: graph
     169              :       INTEGER, INTENT(in)                                :: maxdegree
     170              :       INTEGER, INTENT(out)                               :: ncolors
     171              : 
     172              :       CHARACTER(len=*), PARAMETER :: routineN = 'color_graph_greedy'
     173              : 
     174              :       INTEGER                                            :: color, handle, i, j, newcolor, &
     175              :                                                             nneighbors, nnodes
     176            4 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: color_present
     177              : 
     178            4 :       CALL timeset(routineN, handle)
     179              : 
     180            4 :       ncolors = 0
     181              : 
     182            4 :       nnodes = SIZE(graph)
     183              : 
     184           12 :       ALLOCATE (color_present(maxdegree + 1))
     185              : 
     186           12 :       DO i = 1, nnodes
     187           16 :          color_present(:) = .FALSE.
     188            8 :          IF (ASSOCIATED(graph(i)%vertex%neighbors)) THEN
     189            0 :             nneighbors = SIZE(graph(i)%vertex%neighbors)
     190            0 :             DO j = 1, nneighbors
     191            0 :                color = graph(i)%vertex%neighbors(j)%vertex%color
     192            0 :                IF (color .NE. 0) color_present(color) = .TRUE.
     193              :             END DO
     194              :          END IF
     195            8 :          DO j = 1, maxdegree + 1 !nnodes
     196            8 :             IF (color_present(j) .EQV. .FALSE.) THEN
     197              :                newcolor = j
     198              :                EXIT
     199              :             END IF
     200              :          END DO
     201            8 :          IF (newcolor .GT. ncolors) ncolors = newcolor
     202           12 :          graph(i)%vertex%color = newcolor ! smallest possible
     203              :       END DO
     204              : 
     205            4 :       DEALLOCATE (color_present)
     206              : 
     207            4 :       CALL timestop(handle)
     208              : 
     209            4 :    END SUBROUTINE
     210              : 
     211              :    ! prints the subset info to the screen - useful for vmd visualization
     212              :    ! note that the index starts with '0' and not with '1'
     213              : ! **************************************************************************************************
     214              : !> \brief ...
     215              : !> \param graph ...
     216              : !> \param ncolors ...
     217              : !> \param unit_nr ...
     218              : ! **************************************************************************************************
     219            0 :    SUBROUTINE print_subsets(graph, ncolors, unit_nr)
     220              :       TYPE(vertex_p_type), DIMENSION(:), POINTER         :: graph
     221              :       INTEGER, INTENT(IN)                                :: ncolors, unit_nr
     222              : 
     223              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'print_subsets'
     224              : 
     225              :       INTEGER                                            :: counter, handle, i, j, nnodes
     226              : 
     227            0 :       CALL timeset(routineN, handle)
     228              : 
     229            0 :       IF (unit_nr > 0) THEN
     230              : 
     231            0 :          WRITE (unit_nr, '(T2,A,T10,A)') "Color |", "Molecules in this subset (IDs start from 0)"
     232              : 
     233            0 :          nnodes = SIZE(graph)
     234              : 
     235            0 :          DO i = 1, ncolors
     236            0 :             WRITE (unit_nr, '(T2,I5,1X,A)', ADVANCE='NO') i, "|"
     237            0 :             counter = 0
     238            0 :             DO j = 1, nnodes
     239            0 :                IF (graph(j)%vertex%color .EQ. i) THEN
     240            0 :                   counter = counter + 1
     241            0 :                   IF (MOD(counter, 13) .EQ. 0) THEN
     242            0 :                      counter = counter + 1
     243            0 :                      WRITE (unit_nr, '()') ! line break
     244            0 :                      WRITE (unit_nr, '(6X,A)', ADVANCE='NO') " |" ! indent next line
     245              :                   END IF
     246            0 :                   WRITE (unit_nr, '(I5,1X)', ADVANCE='NO') graph(j)%vertex%id - 1
     247              :                END IF
     248              :             END DO
     249            0 :             WRITE (unit_nr, '()')
     250              :          END DO
     251              : 
     252              :       END IF
     253              : 
     254            0 :       CALL timestop(handle)
     255              : 
     256            0 :    END SUBROUTINE
     257              : 
     258              : ! **************************************************************************************************
     259              : !> \brief ...
     260              : !> \param dsat ...
     261              : !> \param degree ...
     262              : !> \return ...
     263              : ! **************************************************************************************************
     264          136 :    ELEMENTAL FUNCTION kg_get_value(dsat, degree) RESULT(value)
     265              :       INTEGER, INTENT(IN)                                :: dsat, degree
     266              :       INTEGER(KIND=valt)                                 :: value
     267              : 
     268              :       INTEGER, PARAMETER                                 :: huge_4 = 2_int_4**30
     269              : 
     270              : !   INTEGER, PARAMETER                       :: huge_4 = HUGE(0_int_4) ! PR67219 workaround
     271              : ! we actually need a max_heap
     272              : 
     273          136 :       value = (huge_4 - INT(dsat, KIND=int_8))*huge_4 + huge_4 - degree
     274              : 
     275          136 :    END FUNCTION
     276              : 
     277              : ! **************************************************************************************************
     278              : !> \brief ...
     279              : !> \param heap ...
     280              : !> \param graph ...
     281              : ! **************************************************************************************************
     282           37 :    SUBROUTINE kg_cp_heap_fill(heap, graph)
     283              :       TYPE(cp_heap_type), INTENT(INOUT)                  :: heap
     284              :       TYPE(vertex_p_type), DIMENSION(:), INTENT(IN), &
     285              :          POINTER                                         :: graph
     286              : 
     287              :       INTEGER                                            :: i, nnodes
     288           37 :       INTEGER(kind=valt), ALLOCATABLE, DIMENSION(:)      :: values
     289              : 
     290           37 :       nnodes = SIZE(graph)
     291              : 
     292          111 :       ALLOCATE (values(nnodes))
     293              : 
     294          122 :       DO i = 1, nnodes
     295          122 :          values(i) = kg_get_value(0, graph(i)%vertex%degree)
     296              :       END DO
     297              : 
     298              :       ! fill the heap
     299           37 :       CALL cp_heap_fill(heap, values)
     300              : 
     301           37 :       DEALLOCATE (values)
     302              : 
     303           37 :    END SUBROUTINE
     304              : 
     305              : ! **************************************************************************************************
     306              : !> \brief ...
     307              : !> \param heap ...
     308              : !> \param node ...
     309              : ! **************************************************************************************************
     310           51 :    SUBROUTINE kg_update_node(heap, node)
     311              :       TYPE(cp_heap_type)                                 :: heap
     312              :       TYPE(vertex), INTENT(IN), POINTER                  :: node
     313              : 
     314              :       INTEGER                                            :: degree, dsat, id
     315              :       INTEGER(KIND=valt)                                 :: value
     316              : 
     317           51 :       id = node%id
     318              : 
     319              :       ! only update node if not yet colored
     320           51 :       IF (heap%index(id) > 0) THEN
     321              : 
     322           51 :          degree = node%degree
     323           51 :          dsat = node%dsat
     324              : 
     325           51 :          value = kg_get_value(dsat, degree)
     326              : 
     327           51 :          CALL cp_heap_reset_node(heap, id, value)
     328              : 
     329              :       END IF
     330              : 
     331           51 :    END SUBROUTINE
     332              : 
     333              :    ! Subroutine revised by Samuel Andermatt (2013.11)
     334              : ! **************************************************************************************************
     335              : !> \brief ...
     336              : !> \param kg_env ...
     337              : !> \param graph ...
     338              : !> \param ncolors ...
     339              : ! **************************************************************************************************
     340           37 :    SUBROUTINE kg_dsatur(kg_env, graph, ncolors)
     341              :       TYPE(kg_environment_type), POINTER                 :: kg_env
     342              :       TYPE(vertex_p_type), DIMENSION(:), POINTER         :: graph
     343              :       INTEGER(KIND=int_4), INTENT(OUT)                   :: ncolors
     344              : 
     345              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'kg_dsatur'
     346              : 
     347              :       INTEGER                                            :: color_limit, handle, i, j, key, &
     348              :                                                             maxdegree, nneighbors, nnodes
     349              :       INTEGER(KIND=valt)                                 :: value
     350              :       LOGICAL                                            :: found
     351           37 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: color_present
     352              :       TYPE(cp_heap_type)                                 :: heap
     353              :       TYPE(vertex), POINTER                              :: neighbor, this
     354              : 
     355           37 :       CALL timeset(routineN, handle)
     356              : 
     357           37 :       CALL cite_reference(Brelaz1979)
     358              : 
     359           37 :       ncolors = 0
     360           37 :       color_limit = 100
     361           37 :       maxdegree = kg_env%maxdegree
     362           37 :       nnodes = SIZE(graph)
     363              : 
     364           37 :       IF (kg_env%maxdegree .EQ. 0) THEN
     365              :          ! all nodes are disconnected
     366              : 
     367            0 :          ncolors = 1
     368              : 
     369            0 :          DO i = 1, nnodes
     370              :             ! only one color needed to color the graph
     371            0 :             graph(i)%vertex%color = 1
     372              :          END DO
     373              : 
     374              :       ELSE
     375              : 
     376              :          ! allocate and fill the heap
     377           37 :          CALL cp_heap_new(heap, nnodes)
     378              : 
     379           37 :          CALL kg_cp_heap_fill(heap, graph)
     380              : 
     381          122 :          DO WHILE (heap%n > 0)
     382              : 
     383           85 :             CALL cp_heap_pop(heap, key, value, found)
     384              : 
     385           85 :             this => graph(key)%vertex
     386              : 
     387           85 :             nneighbors = 0
     388              : 
     389           85 :             IF (ASSOCIATED(this%neighbors)) THEN
     390           81 :                nneighbors = SIZE(this%neighbors)
     391              :             ELSE !node is isolated
     392            4 :                this%color = 1
     393            4 :                CYCLE
     394              :             END IF
     395          132 :             DO i = 1, this%degree + 1
     396          132 :                IF (this%color_present(i) .EQV. .FALSE.) THEN
     397           81 :                   this%color = i ! smallest possible
     398           81 :                   EXIT
     399              :                END IF
     400              :             END DO
     401           81 :             IF (this%color .GT. ncolors) ncolors = this%color
     402              : 
     403              :             ! update all neighboring nodes
     404          183 :             DO i = 1, nneighbors
     405          102 :                neighbor => this%neighbors(i)%vertex
     406          102 :                IF (neighbor%color_present(this%color)) CYCLE
     407          102 :                neighbor%color_present(this%color) = .TRUE.
     408          102 :                neighbor%dsat = neighbor%dsat + 1
     409          102 :                IF (neighbor%color /= 0) CYCLE
     410          183 :                CALL kg_update_node(heap, neighbor)
     411              : 
     412              :             END DO
     413              : 
     414          118 :             IF (this%color == color_limit) THEN !resize all color_present arrays
     415            0 :                ALLOCATE (color_present(color_limit))
     416            0 :                DO i = 1, nnodes
     417            0 :                   IF (graph(i)%vertex%degree == 0) CYCLE
     418            0 :                   color_present(:) = graph(i)%vertex%color_present
     419            0 :                   DEALLOCATE (graph(i)%vertex%color_present)
     420            0 :                   ALLOCATE (graph(i)%vertex%color_present(color_limit*2))
     421            0 :                   DO j = 1, color_limit
     422            0 :                      graph(i)%vertex%color_present(j) = color_present(j)
     423              :                   END DO
     424            0 :                   DO j = color_limit + 1, 2*color_limit
     425            0 :                      graph(i)%vertex%color_present(j) = .FALSE.
     426              :                   END DO
     427              :                END DO
     428            0 :                DEALLOCATE (color_present)
     429            0 :                color_limit = color_limit*2
     430              :             END IF
     431              : 
     432              :          END DO
     433              : 
     434              :          ! release the heap
     435           37 :          CALL cp_heap_release(heap)
     436              : 
     437              :       END IF
     438              : 
     439           37 :       CALL timestop(handle)
     440              : 
     441           74 :    END SUBROUTINE
     442              : 
     443              : ! **************************************************************************************************
     444              : !> \brief Checks if the color of two nodes can be exchanged legally
     445              : !> \param this ...
     446              : !> \param neighbor ...
     447              : !> \param low_col ...
     448              : !> \param switchable ...
     449              : !> \param ncolors ...
     450              : !> \param color_present ...
     451              : !> \par History
     452              : !>       2013.11 created [Samuel Andermatt]
     453              : !> \author Samuel Andermatt
     454              : ! **************************************************************************************************
     455         6996 :    SUBROUTINE kg_check_switch(this, neighbor, low_col, switchable, ncolors, color_present)
     456              : 
     457              :       TYPE(vertex), POINTER                              :: this, neighbor
     458              :       INTEGER, INTENT(OUT)                               :: low_col
     459              :       LOGICAL                                            :: switchable
     460              :       INTEGER                                            :: ncolors
     461              :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: color_present
     462              : 
     463              :       INTEGER                                            :: i
     464              : 
     465         6996 :       switchable = .TRUE.
     466         6996 :       low_col = ncolors + 1
     467              : 
     468        16218 :       DO i = 1, this%degree
     469         9222 :          IF (this%neighbors(i)%vertex%id == neighbor%id) CYCLE
     470         9222 :          IF (this%neighbors(i)%vertex%color == neighbor%color) THEN
     471            0 :             switchable = .FALSE.
     472            0 :             EXIT
     473              :          END IF
     474              :       END DO
     475         6996 :       IF (.NOT. switchable) RETURN
     476              : 
     477        23214 :       color_present(:) = .FALSE.
     478         6996 :       color_present(neighbor%color) = .TRUE.
     479        16218 :       DO i = 1, neighbor%degree
     480         9222 :          IF (neighbor%neighbors(i)%vertex%id == this%id) CYCLE
     481        16218 :          color_present(neighbor%neighbors(i)%vertex%color) = .TRUE.
     482              :       END DO
     483        16218 :       DO i = 1, this%color
     484        16218 :          IF (.NOT. color_present(i)) THEN
     485         6996 :             low_col = i
     486         6996 :             EXIT
     487              :          END IF
     488              :       END DO
     489              : 
     490              :    END SUBROUTINE
     491              : 
     492              : ! **************************************************************************************************
     493              : !> \brief An algorithm that works on an already colored graph and tries to
     494              : !>          reduce the amount of colors by switching the colors of
     495              : !>          neighboring vertices
     496              : !> \param graph ...
     497              : !> \param ncolors ...
     498              : !> \par History
     499              : !>       2013.11 created [Samuel Andermatt]
     500              : !> \author Samuel Andermatt
     501              : ! **************************************************************************************************
     502           41 :    SUBROUTINE kg_pair_switching(graph, ncolors)
     503              :       TYPE(vertex_p_type), DIMENSION(:), POINTER         :: graph
     504              :       INTEGER                                            :: ncolors
     505              : 
     506              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'kg_pair_switching'
     507              : 
     508              :       INTEGER                                            :: depth, handle, i, iter, j, low_col, &
     509              :                                                             low_col_neigh, maxdepth, &
     510              :                                                             maxiterations, nnodes, partner
     511              :       LOGICAL                                            :: switchable
     512              :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: color_present
     513              :       TYPE(vertex), POINTER                              :: neighbor, this
     514              : 
     515           41 :       CALL timeset(routineN, handle)
     516           41 :       nnodes = SIZE(graph)
     517          123 :       ALLOCATE (color_present(ncolors))
     518              : 
     519              :       !Some default values that should work well
     520           41 :       maxdepth = 4
     521           41 :       maxiterations = 20
     522              : 
     523          861 :       DO iter = 1, maxiterations
     524         4141 :          DO depth = 1, maxdepth !First the nodes with larges colors are attempted to be switched and reduced
     525              :             !Go trough all nodes and try if you can reduce their color by switching with a neighbour
     526        10720 :             DO i = 1, nnodes
     527         7440 :                this => graph(i)%vertex
     528         7440 :                IF (.NOT. ASSOCIATED(this%neighbors)) CYCLE
     529         6480 :                IF (graph(i)%vertex%color < ncolors - depth + 1) CYCLE !Node already has a low color
     530              : 
     531         6163 :                partner = 0
     532         6163 :                low_col = this%color + 1
     533              : 
     534        13719 :                DO j = 1, this%degree
     535         7556 :                   neighbor => this%neighbors(j)%vertex
     536         7556 :                   IF (neighbor%color > this%color) CYCLE
     537         6996 :                   CALL kg_check_switch(this, neighbor, low_col_neigh, switchable, ncolors, color_present)
     538        13159 :                   IF (switchable .AND. low_col_neigh < low_col) THEN
     539         5883 :                      partner = j
     540         5883 :                      low_col = low_col_neigh
     541              :                   END IF
     542              :                END DO
     543         6163 :                IF (partner == 0) CYCLE !Cannot switch favourably with anybody
     544              :                !Switch the nodes
     545         5883 :                this%color = this%neighbors(partner)%vertex%color
     546        10720 :                this%neighbors(partner)%vertex%color = low_col
     547              :             END DO
     548              : 
     549         3280 :             ncolors = 0
     550        11540 :             DO j = 1, nnodes
     551        10720 :                IF (graph(j)%vertex%color > ncolors) THEN
     552         3280 :                   ncolors = graph(j)%vertex%color
     553              :                END IF
     554              :             END DO
     555              :          END DO
     556              :       END DO
     557              : 
     558           41 :       DEALLOCATE (color_present)
     559           41 :       CALL timestop(handle)
     560              : 
     561           41 :    END SUBROUTINE
     562              : 
     563              : ! **************************************************************************************************
     564              : !> \brief ...
     565              : !> \param graph ...
     566              : !> \param valid ...
     567              : ! **************************************************************************************************
     568           41 :    SUBROUTINE check_coloring(graph, valid)
     569              :       TYPE(vertex_p_type), DIMENSION(:), INTENT(in), &
     570              :          POINTER                                         :: graph
     571              :       LOGICAL, INTENT(out)                               :: valid
     572              : 
     573              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'check_coloring'
     574              : 
     575              :       INTEGER                                            :: handle, i, j, nneighbors, nnodes
     576              :       TYPE(vertex), POINTER                              :: neighbor, node
     577              : 
     578           41 :       CALL timeset(routineN, handle)
     579              : 
     580           41 :       valid = .TRUE.
     581           41 :       nnodes = SIZE(graph)
     582              : 
     583          134 :       DO i = 1, nnodes
     584           93 :          node => graph(i)%vertex
     585          134 :          IF (ASSOCIATED(node%neighbors)) THEN
     586           81 :             nneighbors = SIZE(node%neighbors)
     587          183 :             DO j = 1, nneighbors
     588          102 :                neighbor => node%neighbors(j)%vertex
     589          183 :                IF (neighbor%color == node%color) THEN
     590            0 :                   valid = .FALSE.
     591            0 :                   RETURN
     592              :                END IF
     593              :             END DO
     594              :          END IF
     595              :       END DO
     596              : 
     597           41 :       CALL timestop(handle)
     598              : 
     599              :    END SUBROUTINE
     600              : 
     601              : ! **************************************************************************************************
     602              : !> \brief ...
     603              : !> \param graph ...
     604              : ! **************************************************************************************************
     605           41 :    SUBROUTINE deallocate_graph(graph)
     606              :       TYPE(vertex_p_type), DIMENSION(:), POINTER         :: graph
     607              : 
     608              :       INTEGER                                            :: i, nnodes
     609              : 
     610           41 :       nnodes = SIZE(graph)
     611              : 
     612          134 :       DO i = 1, nnodes
     613           93 :          IF (ASSOCIATED(graph(i)%vertex%neighbors)) DEALLOCATE (graph(i)%vertex%neighbors)
     614          134 :          DEALLOCATE (graph(i)%vertex)
     615              :       END DO
     616           41 :       DEALLOCATE (graph)
     617              : 
     618           41 :    END SUBROUTINE
     619              : 
     620              : ! **************************************************************************************************
     621              : !> \brief ...
     622              : !> \param kg_env ...
     623              : !> \param pairs ...
     624              : !> \param ncolors ...
     625              : !> \param color_of_node ...
     626              : ! **************************************************************************************************
     627           41 :    SUBROUTINE kg_vertex_coloring(kg_env, pairs, ncolors, color_of_node)
     628              :       TYPE(kg_environment_type), POINTER                 :: kg_env
     629              :       INTEGER(KIND=int_4), ALLOCATABLE, &
     630              :          DIMENSION(:, :), INTENT(IN)                     :: pairs
     631              :       INTEGER(KIND=int_4), INTENT(OUT)                   :: ncolors
     632              :       INTEGER(KIND=int_4), ALLOCATABLE, DIMENSION(:), &
     633              :          INTENT(out)                                     :: color_of_node
     634              : 
     635              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'kg_vertex_coloring'
     636              : 
     637              :       INTEGER                                            :: color, handle, i, nnodes, unit_nr
     638              :       LOGICAL                                            :: valid
     639              :       TYPE(cp_logger_type), POINTER                      :: logger
     640           41 :       TYPE(vertex_p_type), DIMENSION(:), POINTER         :: graph
     641              : 
     642              : ! get a useful output_unit
     643              : 
     644           82 :       logger => cp_get_default_logger()
     645           41 :       IF (logger%para_env%is_source()) THEN
     646           41 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     647              :       ELSE
     648              :          unit_nr = -1
     649              :       END IF
     650              : 
     651           41 :       CALL timeset(routineN, handle)
     652              : 
     653           41 :       CALL kg_create_graph(kg_env, pairs, graph)
     654              : 
     655           45 :       SELECT CASE (kg_env%coloring_method)
     656              :       CASE (kg_color_greedy)
     657              :          ! simple greedy algorithm
     658           41 :          CALL color_graph_greedy(graph, kg_env%maxdegree, ncolors)
     659              :       CASE (kg_color_dsatur)
     660              :          ! color with max degree of saturation
     661           37 :          CALL kg_dsatur(kg_env, graph, ncolors)
     662              :       CASE DEFAULT
     663           41 :          CPABORT("Coloring method not known.")
     664              :       END SELECT
     665              : 
     666           41 :       CALL kg_pair_switching(graph, ncolors)
     667              : 
     668              :       valid = .FALSE.
     669           41 :       CALL check_coloring(graph, valid)
     670           41 :       IF (.NOT. valid) &
     671            0 :          CPABORT("Coloring not valid.")
     672              : 
     673           41 :       nnodes = SIZE(kg_env%molecule_set)
     674              : 
     675          123 :       ALLOCATE (color_of_node(nnodes))
     676              : 
     677              :       ! gather the subset info in a simple integer array
     678          134 :       DO i = 1, nnodes
     679           93 :          color = graph(i)%vertex%color
     680          134 :          color_of_node(i) = color
     681              :       END DO
     682              : 
     683           41 :       IF (unit_nr > 0) THEN
     684              : 
     685           41 :          WRITE (unit_nr, '(T2,A,A,A)') REPEAT("-", 30), " KG coloring result ", REPEAT("-", 29)
     686              :          IF (.FALSE.) THEN ! more output for VMD
     687              :             WRITE (unit_nr, '()')
     688              :             CALL print_subsets(graph, ncolors, unit_nr)
     689              :             WRITE (unit_nr, '()')
     690              :          END IF
     691           41 :          WRITE (unit_nr, '(T2, A16,50X,I6,1X,A6)') 'Number of colors', ncolors, 'colors'
     692           41 :          IF (valid) WRITE (unit_nr, '(T2, A17,59X,A3)') 'Coloring verified', 'yes'
     693           41 :          WRITE (unit_nr, '(T2,A)') REPEAT("-", 79)
     694              : 
     695              :       END IF
     696              : 
     697           41 :       CALL deallocate_graph(graph)
     698              : 
     699           41 :       CALL timestop(handle)
     700              : 
     701          123 :    END SUBROUTINE
     702              : 
     703            0 : END MODULE
        

Generated by: LCOV version 2.0-1