File: N:\mfix\model\des\pair_manager.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Module: pair_manager                                                !
4     !                                                                      !
5     !  Purpose: maintains a hashtable of pairs of positive 32-bit          !
6     !           integers (meant to represent particle ids).                !
7     !                                                                      !
8     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
9     #include "version.inc"
10     
11     module pair_manager
12     
13       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
14       !                                                                      !
15       !  Type: pair_t                                                        !
16       !                                                                      !
17       !  Purpose: Represents a pair of particle ids.                         !
18       !                                                                      !
19       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
20     
21       type pair_t
22          integer(kind=4) :: ii
23          integer(kind=4) :: jj
24       end type pair_t
25     
26       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
27       !                                                                      !
28       !  Type: hashtable_t                                                   !
29       !                                                                      !
30       !  Purpose: Represents a hashtable of pair_t. Uses open addressing.    !
31       !           Blank values are (0,0); deleted values are (0,1).          !
32       !                                                                      !
33       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
34     
35       type hashtable_t
36          type(pair_t), dimension(:), allocatable :: table
37     
38          ! if table_size/size(table) > 50%, time to resize the hashtable
39          integer :: table_size
40          integer :: current_hash
41       end type hashtable_t
42     
43       public :: init_pairs, reset_pairs, is_pair, get_pair, add_pair, del_pair
44       private :: check_table
45     
46     contains
47     
48       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
49       !                                                                      !
50       !  Subroutine: init_pairs                                              !
51       !                                                                      !
52       !  Purpose: hashtable_t constructor                                    !
53       !                                                                      !
54       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
55     
56       subroutine init_pairs(this)
57         implicit none
58         type(hashtable_t), intent(inout) :: this
59     
60         this%current_hash = 0
61         if (.not.allocated(this%table)) allocate(this%table(0:10006))
62         this%table(:)%ii = 0
63         this%table(:)%jj = 0
64         this%table_size = 0
65     
66       end subroutine init_pairs
67     
68       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
69       !                                                                      !
70       !  Subroutine: reset_pairs                                             !
71       !                                                                      !
72       !  Purpose: Resets the iterator for traversing the hashtable.          !
73       !           Used before calling get_pairs().                           !
74       !                                                                      !
75       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
76     
77       subroutine reset_pairs(this)
78         implicit none
79         type(hashtable_t), intent(inout) :: this
80     
81         this%current_hash = 0
82     
83       end subroutine reset_pairs
84     
85       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
86       !                                                                      !
87       !  Subroutine: check_table                                             !
88       !                                                                      !
89       !  Purpose: Debugging                                                  !
90       !                                                                      !
91       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
92     
93       logical function check_table(this)
94         implicit none
95         type(hashtable_t), intent(in) :: this
96         integer :: nn, blanks, deleted, full
97     
98         if (this%table_size > size(this%table)) then
99            check_table = .false.
100            return
101         endif
102     
103         check_table = .true.
104         ! return
105     
106         blanks = 0
107         deleted = 0
108         full = 0
109         do nn=0, size(this%table)-1
110            if (this%table(nn)%ii > 0 .and. this%table(nn)%jj > 0) then
111               full = full + 1
112            else if (this%table(nn)%ii .eq. 0 .and. this%table(nn)%jj .eq. 0) then
113               blanks = blanks + 1
114            else if (this%table(nn)%ii .eq. 0 .and. this%table(nn)%jj .eq. 1) then
115               deleted = deleted + 1
116            else
117               print *,"SHOULD NEVER OCCUR"
118               check_table = .false.
119               return
120            endif
121         end do
122     
123         if (full .ne. this%table_size) then
124            print *,"SIZE = ",size(this%table)
125            print *,"blanks = ",blanks
126            print *,"deleted = ",deleted
127            print *,"full = ",full
128            print *,"table_size = ",this%table_size
129            check_table = .false.
130         endif
131     
132         if (full+deleted+blanks .ne. size(this%table)) then
133            print *,"SIZE = ",size(this%table)
134            print *,"blanks = ",blanks
135            print *,"deleted = ",deleted
136            print *,"full = ",full
137            print *,"table_size = ",this%table_size
138            check_table = .false.
139         endif
140     
141       end function check_table
142     
143       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
144       !                                                                      !
145       !  Subroutine: get_pair                                                !
146       !                                                                      !
147       !  Purpose: Returns the next pair in the hashtable.  Calling           !
148       !           reset_pairs() starts over again from the beginning.        !
149       !                                                                      !
150       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
151     
152       subroutine get_pair(this,pair)
153         implicit none
154         type(hashtable_t), intent(inout) :: this
155         integer, intent(out) :: pair(2)
156     
157         pair(1) = 0
158         pair(2) = 0
159     
160         do while (this%current_hash < size(this%table))
161            if (0.ne.this%table(this%current_hash)%ii .and. 0.ne.this%table(this%current_hash)%jj) then
162               pair(1) = this%table(this%current_hash)%ii
163               pair(2) = this%table(this%current_hash)%jj
164               this%current_hash = this%current_hash + 1
165               return
166            endif
167            this%current_hash = this%current_hash + 1
168         enddo
169       end subroutine get_pair
170     
171       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
172       !                                                                      !
173       !  Subroutine: is_pair                                                 !
174       !                                                                      !
175       !  Purpose: Returns .true. iff the pair (i0,j0) is in the hashtable.   !
176       !                                                                      !
177       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
178     
179       logical function is_pair(this,i0,j0)
180         implicit none
181         type(hashtable_t), intent(in) :: this
182         integer, intent(in) :: i0, j0
183         integer :: ii, jj, probe_count
184         integer(kind=8) :: hash, init_hash
185     
186         ii = min(i0,j0)
187         jj = max(i0,j0)
188     
189         if (ii < 1 .or. jj < 1) then
190            print *,"invalid pair: ",ii,jj
191            ERROR_STOP __LINE__
192         endif
193     
194         ! assign ii to hash to convert to 64-bit
195         probe_count = 1
196         hash = mod(ii+jj*jj+probe_count*probe_count,size(this%table))
197         if (hash < 0) hash = hash+size(this%table)
198         init_hash = hash
199         ! print *,"INIT HASH IS =",hash," TABLE IS ",this%table_size,"/",size(this%table)
200     
201         do
202            if (this%table(hash)%ii .eq. ii .and. this%table(hash)%jj .eq. jj) then
203               is_pair = .true.
204               return
205            endif
206            if (this%table(hash)%ii .eq. 0 .and. this%table(hash)%jj .eq. 0) then
207               is_pair = .false.
208               return
209            endif
210            probe_count = probe_count + 1
211            hash = mod(hash+probe_count*probe_count,size(this%table))
212            if (hash < 0) hash = hash+size(this%table)
213            if (hash .eq. init_hash) exit
214         enddo
215     
216         print *,"loop in hash addressing, this should not occur"
217         ERROR_STOP __LINE__
218     
219       end function is_pair
220     
221       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
222       !                                                                      !
223       !  Subroutine: add_pair                                                !
224       !                                                                      !
225       !  Purpose: Adds (i0,j0) to the hashtable; resizes if necessary.       !
226       !                                                                      !
227       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
228     
229       recursive subroutine add_pair(this,i0,j0)
230         implicit none
231         type(hashtable_t), intent(inout) :: this
232         integer, intent(in) :: i0,j0
233         integer :: ii, jj, nn, old_size, old_tablesize, probe_count
234         integer(kind=8) :: hash, init_hash
235         type(pair_t), dimension(:), allocatable :: table_tmp
236     
237         if (i0 < 1 .or. j0 < 1) then
238            print *,"invalid pair: ",i0,j0
239            ERROR_STOP __LINE__
240         endif
241     
242         if (size(this%table) < 2*this%table_size ) then
243            old_size = size(this%table)
244            old_tablesize = this%table_size
245            allocate(table_tmp(0:old_size-1))
246            if (size(table_tmp).ne.old_size) then
247               print *,"size = ",size(table_tmp)
248               print *,"old_size = ",old_size
249               ERROR_STOP __LINE__
250            endif
251            table_tmp(0:old_size-1) = this%table(0:old_size-1)
252     
253            deallocate(this%table)
254            allocate(this%table(0:2*old_size))
255            this%table(:)%ii = 0
256            this%table(:)%jj = 0
257            this%table_size = 0
258            do nn=0, old_size-1
259               if ( table_tmp(nn)%ii .ne. 0 .and. table_tmp(nn)%jj .ne. 0) then
260                  call add_pair(this,table_tmp(nn)%ii,table_tmp(nn)%jj)
261               endif
262            enddo
263            if (this%table_size.ne.old_tablesize) then
264               print *,"size = ",this%table_size
265               print *,"old_size = ",old_tablesize
266               ERROR_STOP __LINE__
267            endif
268            deallocate(table_tmp)
269         endif
270     
271         ii = min(i0,j0)
272         jj = max(i0,j0)
273     
274         if (ii < 1 .or. jj < 1) then
275            print *,"invalid pair: ",ii,jj
276            ERROR_STOP __LINE__
277         endif
278     
279         ! assign ii to hash to convert to 64-bit
280         probe_count = 1
281         hash = mod(ii+jj*jj+probe_count*probe_count,size(this%table))
282         if (hash < 0) hash = hash+size(this%table)
283         init_hash = hash
284     
285         do
286            if (this%table(hash)%ii .eq. ii .and. this%table(hash)%jj .eq. jj) then
287               ! already in table
288               return
289            endif
290            if (this%table(hash)%ii .eq. 0 .or. this%table(hash)%jj .eq. 0) then
291               this%table(hash)%ii = ii
292               this%table(hash)%jj = jj
293               this%table_size = this%table_size + 1
294               ! if(.not. check_table(this)) ERROR_STOP __LINE__
295               return
296            endif
297            probe_count = probe_count + 1
298            hash = mod(hash+probe_count*probe_count,size(this%table))
299            if (hash < 0) hash = hash+size(this%table)
300            if (hash .eq. init_hash) exit
301         enddo
302     
303         print *,"loop in hash addressing, this should not occur.  maybe hash table is full"
304         ERROR_STOP __LINE__
305     
306       end subroutine add_pair
307     
308       !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
309       !                                                                      !
310       !  Subroutine: del_pair                                                !
311       !                                                                      !
312       !  Purpose: Removes (i0,j0) from the hashtable (if it exists).         !
313       !                                                                      !
314       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
315     
316       subroutine del_pair(this,i0,j0)
317         implicit none
318         type(hashtable_t), intent(inout) :: this
319         integer, intent(in) :: i0,j0
320         integer :: ii, jj, probe_count
321         integer(kind=8) :: hash, init_hash
322     
323         ii = min(i0,j0)
324         jj = max(i0,j0)
325     
326         if (ii < 1 .or. jj < 1) then
327            print *,"invalid pair: ",ii,jj
328            ERROR_STOP __LINE__
329         endif
330     
331         ! assign ii to hash to convert to 64-bit
332         probe_count = 1
333         hash = mod(ii+jj*jj+probe_count*probe_count,size(this%table))
334         if (hash < 0) hash = hash+size(this%table)
335         init_hash = hash
336     
337         do
338            if (this%table(hash)%ii .eq. 0 .and. this%table(hash)%jj .eq. 0) then
339               ! if(.not. check_table(this)) ERROR_STOP __LINE__
340               return
341            endif
342            if (this%table(hash)%ii .eq. ii .and. this%table(hash)%jj .eq. jj) then
343               ! 0,1 signifies DELETED hash entry
344               this%table(hash)%ii = 0
345               this%table(hash)%jj = 1
346               this%table_size = this%table_size - 1
347               ! if(.not. check_table(this)) ERROR_STOP __LINE__
348               return
349            endif
350            probe_count = probe_count + 1
351            hash = mod(hash+probe_count*probe_count,size(this%table))
352            if (hash < 0) hash = hash+size(this%table)
353            if (hash .eq. init_hash) exit
354         enddo
355     
356         ! if(.not. check_table(this)) ERROR_STOP __LINE__
357     
358         print *,"loop in hash addressing. must be a lot of DELETED entries:  ",this%table_size,"/",size(this%table)
359     
360       end subroutine del_pair
361     
362     end module pair_manager
363