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