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 uses a combination of graphs and hashing to determine if two molecules
10 : !> are topologically equivalent, and if so, finds the one by one mapping
11 : !> \note
12 : !> the graph isomorphism being solved is a computationally hard one
13 : !> and can not be solved in polynomial time in the general case
14 : !> http://mathworld.wolfram.com/IsomorphicGraphs.html
15 : !> the problem arises if many atoms are topologically equivalent
16 : !> the current algorithm is able to solve the problem for benzene (C6H6)
17 : !> but not for a fullerene (C60). Large systems are not really a problem (JAC).
18 : !> as almost all atoms are topologically unique.
19 : !> \par History
20 : !> 09.2006 [Joost VandeVondele]
21 : !> \author Joost VandeVondele
22 : ! **************************************************************************************************
23 : MODULE graphcon
24 :
25 : USE util, ONLY: sort
26 : #include "./base/base_uses.f90"
27 :
28 : IMPLICIT NONE
29 :
30 : PRIVATE
31 : PUBLIC :: vertex, graph_type, reorder_graph, hash_molecule
32 :
33 : ! a molecule is an array of vertices, each vertex has a kind
34 : ! and a list of edges (bonds).
35 : ! (the number is the index of the other vertex in the array that builds the molecule)
36 : ! **************************************************************************************************
37 : TYPE graph_type
38 : TYPE(vertex), POINTER, DIMENSION(:) :: graph
39 : END TYPE graph_type
40 :
41 : ! **************************************************************************************************
42 : TYPE vertex
43 : INTEGER :: kind
44 : INTEGER, POINTER, DIMENSION(:) :: bonds
45 : END TYPE vertex
46 :
47 : ! **************************************************************************************************
48 : TYPE class
49 : INTEGER, DIMENSION(:), POINTER :: reference
50 : INTEGER, DIMENSION(:), POINTER :: unordered
51 : INTEGER :: kind
52 : INTEGER :: Nele
53 : LOGICAL :: first
54 : INTEGER, DIMENSION(:), POINTER :: order
55 : INTEGER, DIMENSION(:), POINTER :: q
56 : END TYPE class
57 :
58 : ! **************************************************************************************************
59 : TYPE superclass
60 : INTEGER :: Nele
61 : INTEGER, DIMENSION(:), POINTER :: classes
62 : END TYPE
63 :
64 : CONTAINS
65 :
66 : ! **************************************************************************************************
67 : !> \brief hashes a molecule to a number. Molecules that are the (topologically) the same
68 : !> have the same hash. However, there is a small chance that molecules with the same hash
69 : !> are different
70 : !> \param reference IN : molecule with atomic kinds and bonds
71 : !> \param kind_ref OUT : an atomic hash which is the same for topologically equivalent atoms
72 : !> \param hash OUT : a hash which is the same for topologically equivalent molecules
73 : !> \par History
74 : !> 09.2006 created [Joost VandeVondele]
75 : !> \note
76 : !> Although relatively fast in general, might be quadratic with molecule size for
77 : !> some systems (e.g. linear alkanes)
78 : ! **************************************************************************************************
79 1468 : SUBROUTINE hash_molecule(reference, kind_ref, hash)
80 : TYPE(vertex), DIMENSION(:), INTENT(IN) :: reference
81 : INTEGER, DIMENSION(:), INTENT(OUT) :: kind_ref
82 : INTEGER, INTENT(OUT) :: hash
83 :
84 : INTEGER :: I, Ihash, N, Nclasses, Nclasses_old, &
85 : old_class
86 1468 : INTEGER, ALLOCATABLE, DIMENSION(:) :: index, kind_new
87 :
88 1468 : N = SIZE(kind_ref)
89 5872 : ALLOCATE (kind_new(N), INDEX(N))
90 10174 : kind_ref = reference%kind
91 : Nclasses_old = 0
92 3230 : DO Ihash = 1, N
93 : ! generate a hash based on the the kind of each atom and the kind of its bonded atoms
94 23592 : DO I = 1, N
95 100056 : kind_new(I) = hash_kind(kind_ref(I), kind_ref(reference(I)%bonds))
96 : END DO
97 23592 : kind_ref = kind_new
98 : ! find the number of equivalent atoms
99 3150 : CALL sort(kind_new, N, index)
100 3150 : Nclasses = 1
101 3150 : old_class = kind_new(1)
102 20442 : DO i = 2, N
103 20442 : IF (kind_new(I) .NE. old_class) THEN
104 6592 : Nclasses = Nclasses + 1
105 6592 : old_class = kind_new(I)
106 : END IF
107 : END DO
108 : ! if we have not generated new classes, we have presumably found all equivalence classes
109 3150 : IF (Nclasses == Nclasses_old) EXIT
110 3230 : Nclasses_old = Nclasses
111 : ! write(*,*) "Classes",Ihash, Nclasses
112 : END DO
113 : ! hash (sorted) kinds to a molecular hash
114 1468 : hash = joaat_hash_i(kind_new)
115 1468 : DEALLOCATE (kind_new, index)
116 1468 : END SUBROUTINE hash_molecule
117 :
118 : ! **************************************************************************************************
119 : !> \brief If two molecules are topologically the same, finds the ordering that maps
120 : !> the unordered one on the ordered one.
121 : !> \param reference molecular description (see type definition)
122 : !> \param unordered molecular description (see type definition)
123 : !> \param order the mapping reference=order(unordred) if matches=.TRUE.
124 : !> undefined if matches=.FALSE.
125 : !> \param matches .TRUE. = the ordering was found
126 : !> \par History
127 : !> 09.2006 created [Joost VandeVondele]
128 : !> \note
129 : !> See not at the top of the file about why this algorithm might consider
130 : !> molecules with a large number of equivalent atoms as different
131 : !> despite the fact that an ordering could exist for which they are the same
132 : ! **************************************************************************************************
133 1150 : SUBROUTINE reorder_graph(reference, unordered, order, matches)
134 : TYPE(vertex), DIMENSION(:), INTENT(IN) :: reference, unordered
135 : INTEGER, DIMENSION(:), INTENT(OUT) :: order
136 : LOGICAL, INTENT(OUT) :: matches
137 :
138 : INTEGER, PARAMETER :: max_tries = 1000000
139 :
140 : INTEGER :: hash_re, hash_un, I, Iclass, iele, &
141 : isuperclass, itries, J, N, Nclasses, &
142 : Nele, old_class
143 1150 : INTEGER, ALLOCATABLE, DIMENSION(:) :: class_of_atom, index_ref, index_un, &
144 1150 : kind_ref, kind_ref_ordered, kind_un, &
145 1150 : kind_un_ordered, superclass_of_atom
146 1150 : TYPE(class), ALLOCATABLE, DIMENSION(:) :: classes
147 1150 : TYPE(superclass), ALLOCATABLE, DIMENSION(:) :: superclasses
148 :
149 : ! allows for worst case matching of two benzenes ... (6!)*(6!)/6=86400
150 : ! with some margin for other molecules
151 : ! molecules with no symmetry e.g. JAC need less than 500 tries
152 : ! catch the cases where the molecules are trivially different
153 :
154 1150 : IF (SIZE(reference) .NE. SIZE(unordered)) THEN
155 38 : matches = .FALSE.
156 : RETURN
157 : END IF
158 :
159 : ! catch the case where the molecules are already in the right order
160 1112 : N = SIZE(order)
161 7644 : order = (/(i, i=1, N)/)
162 1112 : IF (matrix_equal(reference, unordered, order)) THEN
163 1088 : matches = .TRUE.
164 1088 : RETURN
165 : END IF
166 :
167 : ! determine the kind of each atom, and the hash of the whole molecule
168 : ALLOCATE (kind_ref(N), kind_un(N), index_ref(N), index_un(N), &
169 : kind_ref_ordered(N), kind_un_ordered(N), &
170 240 : class_of_atom(N), superclass_of_atom(N))
171 24 : CALL hash_molecule(reference, kind_ref, hash_re)
172 24 : CALL hash_molecule(unordered, kind_un, hash_un)
173 24 : IF (hash_re .NE. hash_un) THEN
174 18 : matches = .FALSE.
175 18 : RETURN
176 : END IF
177 :
178 : ! generate the classes of equivalent atoms, i.e. the groups of atoms of the same topological kind
179 62 : kind_ref_ordered(:) = kind_ref
180 6 : CALL sort(kind_ref_ordered, N, index_ref)
181 62 : kind_un_ordered(:) = kind_un
182 6 : CALL sort(kind_un_ordered, N, index_un)
183 62 : IF (ANY(kind_ref_ordered .NE. kind_un_ordered)) THEN
184 0 : matches = .FALSE.
185 0 : RETURN
186 : END IF
187 :
188 : ! count different classes, assign their kinds, and the number of elements
189 6 : Nclasses = 1
190 6 : old_class = kind_ref_ordered(1)
191 56 : DO i = 2, N
192 56 : IF (kind_ref_ordered(I) .NE. old_class) THEN
193 20 : Nclasses = Nclasses + 1
194 20 : old_class = kind_ref_ordered(I)
195 : END IF
196 : END DO
197 18 : ALLOCATE (classes(Nclasses))
198 6 : classes(1)%kind = kind_ref_ordered(1)
199 6 : Nclasses = 1
200 6 : classes(1)%Nele = 1
201 56 : DO i = 2, N
202 56 : IF (kind_ref_ordered(I) .NE. classes(Nclasses)%kind) THEN
203 20 : Nclasses = Nclasses + 1
204 20 : classes(Nclasses)%kind = kind_ref_ordered(I)
205 20 : classes(Nclasses)%Nele = 1
206 : ELSE
207 30 : classes(Nclasses)%Nele = classes(Nclasses)%Nele + 1
208 : END IF
209 : END DO
210 :
211 : ! assign the atoms to their classes
212 6 : iele = 0
213 32 : DO I = 1, Nclasses
214 26 : Nele = classes(I)%Nele
215 78 : ALLOCATE (classes(I)%reference(Nele))
216 52 : ALLOCATE (classes(I)%unordered(Nele))
217 82 : DO J = 1, Nele
218 56 : iele = iele + 1
219 56 : classes(I)%reference(J) = index_ref(iele)
220 82 : classes(I)%unordered(J) = index_un(iele)
221 : END DO
222 138 : class_of_atom(classes(I)%reference) = I
223 52 : ALLOCATE (classes(I)%order(Nele))
224 52 : ALLOCATE (classes(I)%q(Nele))
225 138 : classes(I)%order = (/(J, J=1, Nele)/)
226 32 : classes(I)%first = .TRUE.
227 : END DO
228 :
229 : ! find which groups of classes (superclasses) that can be solved independently.
230 : ! only classes with more than one element that are connected need to be reordered simultaniously
231 :
232 : ! find these connected components in a recursive way
233 62 : superclass_of_atom = -1
234 6 : isuperclass = 0
235 62 : DO I = 1, N
236 : ! this atom belongs to a class with several equivalent atoms, and has not yet been found
237 62 : IF (superclass_of_atom(I) .EQ. -1 .AND. classes(class_of_atom(I))%Nele > 1) THEN
238 6 : isuperclass = isuperclass + 1
239 6 : CALL spread_superclass(I, isuperclass, superclass_of_atom, class_of_atom, classes, reference)
240 : END IF
241 : END DO
242 :
243 : ! put classes into superclasses
244 18 : ALLOCATE (superclasses(isuperclass))
245 12 : superclasses%Nele = 0
246 32 : DO I = 1, Nclasses
247 26 : J = superclass_of_atom(classes(I)%reference(1))
248 32 : IF (J > 0) superclasses(J)%Nele = superclasses(J)%Nele + 1
249 : END DO
250 12 : DO I = 1, isuperclass
251 18 : ALLOCATE (superclasses(I)%classes(superclasses(I)%Nele))
252 12 : superclasses(I)%Nele = 0
253 : END DO
254 32 : DO I = 1, Nclasses
255 26 : J = superclass_of_atom(classes(I)%reference(1))
256 32 : IF (J > 0) THEN
257 14 : superclasses(J)%Nele = superclasses(J)%Nele + 1
258 14 : superclasses(J)%classes(superclasses(J)%Nele) = I
259 : END IF
260 : END DO
261 :
262 : ! write(*,*) "Class generation time",t2-t1
263 : ! WRITE(*,*) "Nclasses, max size, total-non-1 ",Nclasses,MAXVAL(classes%Nele),COUNT(classes%Nele>1)
264 : ! write(*,*) "isuperclass ",isuperclass
265 :
266 : ! assign the order array to their initial value
267 32 : DO Iclass = 1, Nclasses
268 200 : order(classes(Iclass)%unordered) = classes(Iclass)%reference(classes(Iclass)%order)
269 : END DO
270 :
271 : ! reorder the atoms superclass after superclass
272 6 : itries = 0
273 12 : DO I = 1, isuperclass
274 : DO
275 87434 : itries = itries + 1
276 :
277 : ! assign the current order
278 262324 : DO iclass = 1, superclasses(I)%Nele
279 174890 : J = superclasses(I)%classes(iclass)
280 3409744 : order(classes(J)%unordered) = classes(J)%reference(classes(J)%order)
281 : END DO
282 :
283 : ! check for matches within this superclass only, be happy if we have a match
284 87434 : matches = matrix_superclass_equal(reference, unordered, order, superclasses(I), classes)
285 87434 : IF (itries > max_tries) THEN
286 0 : WRITE (*, *) "Could not find the 1-to-1 mapping to prove graph isomorphism"
287 0 : WRITE (*, *) "Reordering failed, assuming these molecules are different"
288 0 : EXIT
289 : END IF
290 87434 : IF (matches) EXIT
291 :
292 : ! generate next permutation within this superclass
293 87554 : DO iclass = 1, superclasses(I)%Nele
294 87554 : J = superclasses(I)%classes(iclass)
295 : CALL all_permutations(classes(J)%order, classes(J)%Nele, &
296 87554 : classes(J)%q, classes(J)%first)
297 87554 : IF (.NOT. classes(J)%first) EXIT
298 : END DO
299 :
300 : ! we are back at the original permutation so we're unable to match this superclass.
301 87428 : IF (iclass .EQ. superclasses(I)%Nele .AND. &
302 6 : classes(superclasses(I)%classes(superclasses(I)%Nele))%first) EXIT
303 : END DO
304 : ! failed in this superblock, can exit now
305 12 : IF (.NOT. matches) EXIT
306 : END DO
307 :
308 : ! the final check, just to be sure
309 6 : matches = matrix_equal(reference, unordered, order)
310 :
311 32 : DO Iclass = 1, Nclasses
312 26 : DEALLOCATE (classes(Iclass)%reference)
313 26 : DEALLOCATE (classes(Iclass)%unordered)
314 26 : DEALLOCATE (classes(Iclass)%order)
315 32 : DEALLOCATE (classes(Iclass)%q)
316 : END DO
317 6 : DEALLOCATE (classes)
318 12 : DO I = 1, isuperclass
319 12 : DEALLOCATE (superclasses(I)%classes)
320 : END DO
321 6 : DEALLOCATE (superclasses)
322 1150 : END SUBROUTINE reorder_graph
323 :
324 : ! **************************************************************************************************
325 : !> \brief spreads the superclass over all atoms of this class and all their bonded atoms
326 : !> provided that the latter belong to a class which contains more than one element
327 : !> \param I ...
328 : !> \param isuperclass ...
329 : !> \param superclass_of_atom ...
330 : !> \param class_of_atom ...
331 : !> \param classes ...
332 : !> \param reference ...
333 : !> \par History
334 : !> 09.2006 created [Joost VandeVondele]
335 : ! **************************************************************************************************
336 274 : RECURSIVE SUBROUTINE spread_superclass(I, isuperclass, superclass_of_atom, class_of_atom, &
337 274 : classes, reference)
338 : INTEGER, INTENT(IN) :: i, isuperclass
339 : INTEGER, DIMENSION(:), INTENT(INOUT) :: superclass_of_atom
340 : INTEGER, DIMENSION(:), INTENT(IN) :: class_of_atom
341 : TYPE(class), DIMENSION(:), INTENT(IN) :: classes
342 : TYPE(vertex), DIMENSION(:), INTENT(IN) :: reference
343 :
344 : INTEGER :: J
345 :
346 274 : IF (superclass_of_atom(I) .EQ. -1 .AND. classes(class_of_atom(I))%Nele > 1) THEN
347 44 : superclass_of_atom(I) = isuperclass
348 228 : DO J = 1, classes(class_of_atom(I))%Nele
349 : CALL spread_superclass(classes(class_of_atom(I))%reference(J), isuperclass, &
350 228 : superclass_of_atom, class_of_atom, classes, reference)
351 : END DO
352 128 : DO J = 1, SIZE(reference(I)%bonds)
353 : CALL spread_superclass(reference(I)%bonds(J), isuperclass, &
354 128 : superclass_of_atom, class_of_atom, classes, reference)
355 : END DO
356 : END IF
357 274 : END SUBROUTINE spread_superclass
358 :
359 : ! **************************************************************************************************
360 : !> \brief determines of the vertices of this superclass have the same edges
361 : !> \param reference ...
362 : !> \param unordered ...
363 : !> \param order ...
364 : !> \param super ...
365 : !> \param classes ...
366 : !> \return ...
367 : !> \par History
368 : !> 09.2006 created [Joost VandeVondele]
369 : ! **************************************************************************************************
370 87434 : FUNCTION matrix_superclass_equal(reference, unordered, order, super, classes) RESULT(res)
371 : TYPE(vertex), DIMENSION(:), INTENT(IN) :: reference, unordered
372 : INTEGER, DIMENSION(:), INTENT(IN) :: order
373 : TYPE(superclass), INTENT(IN) :: super
374 : TYPE(class), DIMENSION(:), INTENT(IN) :: classes
375 : LOGICAL :: res
376 :
377 : INTEGER :: I, iclass, iele, J
378 :
379 : ! I is the atom in the unordered set
380 :
381 87574 : loop: DO iclass = 1, super%Nele
382 106190 : DO iele = 1, classes(super%classes(iclass))%Nele
383 106044 : I = classes(super%classes(iclass))%unordered(iele)
384 : res = (reference(order(I))%kind == unordered(I)%kind .AND. &
385 106044 : SIZE(reference(order(I))%bonds) == SIZE(unordered(I)%bonds))
386 140 : IF (res) THEN
387 124734 : DO J = 1, SIZE(reference(order(I))%bonds)
388 212506 : IF (ALL(reference(order(I))%bonds(:) .NE. order(unordered(I)%bonds(J)))) THEN
389 : res = .FALSE.
390 : EXIT loop
391 : END IF
392 : END DO
393 : ELSE
394 : EXIT loop
395 : END IF
396 : END DO
397 : END DO loop
398 87434 : END FUNCTION matrix_superclass_equal
399 :
400 : ! **************************************************************************************************
401 : !> \brief determines of the vertices of the full set is equal, i.e.
402 : !> we have the same connectivity graph
403 : !> \param reference ...
404 : !> \param unordered ...
405 : !> \param order ...
406 : !> \return ...
407 : !> \par History
408 : !> 09.2006 created [Joost VandeVondele]
409 : ! **************************************************************************************************
410 1118 : FUNCTION matrix_equal(reference, unordered, order) RESULT(res)
411 : TYPE(vertex), DIMENSION(:), INTENT(IN) :: reference, unordered
412 : INTEGER, DIMENSION(:), INTENT(IN) :: order
413 : LOGICAL :: res
414 :
415 : INTEGER :: I, J
416 :
417 4366 : loop: DO I = 1, SIZE(reference)
418 : res = (reference(order(I))%kind == unordered(I)%kind .AND. &
419 3272 : SIZE(reference(order(I))%bonds) == SIZE(unordered(I)%bonds))
420 1094 : IF (res) THEN
421 7564 : DO J = 1, SIZE(reference(order(I))%bonds)
422 8692 : IF (ALL(reference(order(I))%bonds(:) .NE. order(unordered(I)%bonds(J)))) THEN
423 : res = .FALSE.
424 : EXIT loop
425 : END IF
426 : END DO
427 : ELSE
428 : EXIT loop
429 : END IF
430 : END DO loop
431 1118 : END FUNCTION matrix_equal
432 :
433 : ! **************************************************************************************************
434 : !> \brief creates a hash for an atom based on its own kind and on the kinds
435 : !> of its bonded neighbors
436 : !> \param me ...
437 : !> \param bonds ...
438 : !> \return ...
439 : !> \par History
440 : !> 09.2006 created [Joost VandeVondele]
441 : !> \note
442 : !> bonds are sorted so that the order of neighbors appearing in the bonded list
443 : !> is not important
444 : ! **************************************************************************************************
445 20442 : FUNCTION hash_kind(me, bonds) RESULT(res)
446 : INTEGER, INTENT(IN) :: me
447 : INTEGER, DIMENSION(:), INTENT(IN) :: bonds
448 : INTEGER :: res
449 :
450 : INTEGER :: I, N
451 20442 : INTEGER, ALLOCATABLE, DIMENSION(:) :: index, ordered_bonds
452 :
453 20442 : N = SIZE(bonds)
454 102130 : ALLOCATE (ordered_bonds(N + 1), INDEX(N))
455 58674 : DO I = 1, N
456 58674 : ordered_bonds(I) = bonds(I)
457 : END DO
458 20442 : ordered_bonds(N + 1) = me
459 : ! N: only sort the bonds, not me
460 20442 : CALL sort(ordered_bonds, N, index)
461 20442 : res = joaat_hash_i(ordered_bonds)
462 20442 : END FUNCTION hash_kind
463 :
464 : ! **************************************************************************************************
465 : !> \brief generates the hash of an array of integers and the index in the table
466 : !> \param key an integer array of any length
467 : !> \return ...
468 : !> \par History
469 : !> 09.2006 created [Joost VandeVondele]
470 : !> \note
471 : !> http://en.wikipedia.org/wiki/Hash_table
472 : !> http://www.burtleburtle.net/bob/hash/doobs.html
473 : !> However, since fortran doesn't have an unsigned 4 byte int
474 : !> we compute it using an integer with the appropriate range
475 : !> we return already the index in the table as a final result
476 : ! **************************************************************************************************
477 21910 : FUNCTION joaat_hash_i(key) RESULT(hash_index)
478 : INTEGER, DIMENSION(:), INTENT(IN) :: key
479 : INTEGER :: hash_index
480 :
481 : INTEGER, PARAMETER :: k64 = SELECTED_INT_KIND(10)
482 : INTEGER(KIND=k64), PARAMETER :: b32 = 2_k64**32 - 1_k64
483 :
484 : INTEGER :: i
485 : INTEGER(KIND=k64) :: hash
486 :
487 21910 : hash = 0_k64
488 89290 : DO i = 1, SIZE(key)
489 67380 : hash = IAND(hash + IBITS(key(i), 0, 8), b32)
490 67380 : hash = IAND(hash + IAND(ISHFT(hash, 10), b32), b32)
491 67380 : hash = IAND(IEOR(hash, IAND(ISHFT(hash, -6), b32)), b32)
492 67380 : hash = IAND(hash + IBITS(key(i), 8, 8), b32)
493 67380 : hash = IAND(hash + IAND(ISHFT(hash, 10), b32), b32)
494 67380 : hash = IAND(IEOR(hash, IAND(ISHFT(hash, -6), b32)), b32)
495 67380 : hash = IAND(hash + IBITS(key(i), 16, 8), b32)
496 67380 : hash = IAND(hash + IAND(ISHFT(hash, 10), b32), b32)
497 67380 : hash = IAND(IEOR(hash, IAND(ISHFT(hash, -6), b32)), b32)
498 67380 : hash = IAND(hash + IBITS(key(i), 24, 8), b32)
499 67380 : hash = IAND(hash + IAND(ISHFT(hash, 10), b32), b32)
500 89290 : hash = IAND(IEOR(hash, IAND(ISHFT(hash, -6), b32)), b32)
501 : END DO
502 21910 : hash = IAND(hash + IAND(ISHFT(hash, 3), b32), b32)
503 21910 : hash = IAND(IEOR(hash, IAND(ISHFT(hash, -11), b32)), b32)
504 21910 : hash = IAND(hash + IAND(ISHFT(hash, 15), b32), b32)
505 : ! hash is the real 32bit hash value of the string,
506 : ! hash_index is an index in the hash_table
507 21910 : hash_index = INT(MOD(hash, INT(HUGE(hash_index), KIND=k64)), KIND=KIND(hash_index))
508 21910 : END FUNCTION joaat_hash_i
509 :
510 : !===ACM Algorithm 323, Generation of Permutations in Lexicographic
511 : ! Order (G6) by R. J. Ord-Smith, CACM 11 (Feb. 1968):117
512 : ! Original Algorithm modified via Certification by I.M. Leitch,
513 : ! 17 March 1969.
514 : ! Algol to Fortran 77 by H.D.Knoble <hdkLESS at SPAM psu dot edu>,
515 : ! May 1995.
516 : ! x = initial values (/1...n/), first=.TRUE.
517 : ! q = scratch
518 : ! first = .TRUE. if you're back at the original order
519 : ! **************************************************************************************************
520 : !> \brief ...
521 : !> \param x ...
522 : !> \param n ...
523 : !> \param q ...
524 : !> \param first ...
525 : ! **************************************************************************************************
526 87554 : SUBROUTINE all_permutations(x, n, q, first)
527 : INTEGER :: n, x(n), q(n)
528 : LOGICAL :: first
529 :
530 : INTEGER :: k, m, t
531 :
532 87554 : IF (n == 1) RETURN
533 87554 : IF (first) THEN
534 134 : first = .FALSE.
535 764 : DO m = 1, n - 1
536 764 : q(m) = n
537 : END DO
538 : END IF
539 87554 : IF (q(n - 1) .EQ. n) THEN
540 43780 : q(n - 1) = n - 1
541 43780 : t = x(n)
542 43780 : x(n) = x(n - 1)
543 43780 : x(n - 1) = t
544 43780 : RETURN
545 : END IF
546 106630 : DO k = n - 1, 1, -1
547 106630 : IF (q(k) .EQ. k) THEN
548 62856 : q(k) = n
549 : ELSE
550 : go to 1
551 : END IF
552 : END DO
553 126 : first = .TRUE.
554 126 : k = 1
555 43774 : GOTO 2
556 43648 : 1 m = q(k)
557 43648 : t = x(m)
558 43648 : x(m) = x(k)
559 43648 : x(k) = t
560 43648 : q(k) = m - 1
561 43648 : k = k + 1
562 43774 : 2 m = n
563 43774 : t = x(m)
564 43774 : x(m) = x(k)
565 43774 : x(k) = t
566 43774 : m = m - 1
567 43774 : k = k + 1
568 47540 : DO WHILE (k .LT. m)
569 3766 : t = x(m)
570 3766 : x(m) = x(k)
571 3766 : x(k) = t
572 3766 : m = m - 1
573 3766 : k = k + 1
574 : END DO
575 : END SUBROUTINE
576 0 : END MODULE graphcon
577 :
|
