Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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
46 : INTEGER :: color
47 : INTEGER :: degree ! degree (number of neighbors)
48 : INTEGER :: dsat ! 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
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 123 : 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 132 : ALLOCATE (graph(imol)%vertex%neighbors(nneighbors))
114 44 : graph(imol)%vertex%degree = nneighbors
115 : 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 111 : ALLOCATE (graph(imol)%vertex%neighbors(nneighbors))
125 37 : graph(imol)%vertex%degree = nneighbors
126 : 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 4 : 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
|
