223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

added pInts and pReals

This commit is contained in:
Martin Diehl 2012-02-09 15:55:04 +00:00
parent 3228cf563c
commit 724ec040a2
1 changed files with 37 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
code/kdtree2.f90
View File

@ -20,8 +20,8 @@ module kdtree2_priority_queue_module
! !
type kdtree2_result type kdtree2_result
! a pair of distances, indexes ! a pair of distances, indexes
real(pReal) :: dis !=0.0 real(pReal) :: dis = 0.0_pReal
integer(pInt) :: idx !=-1 Initializers cause some bugs in compilers. integer(pInt) :: idx = -1_pInt !Initializers cause some bugs in compilers.
end type kdtree2_result end type kdtree2_result
! !
! A heap-based priority queue lets one efficiently implement the following ! A heap-based priority queue lets one efficiently implement the following
@ -267,8 +267,8 @@ bigloop: do
! move last element to first ! move last element to first
! !
a%elems(1) = a%elems(a%heap_size) a%elems(1) = a%elems(a%heap_size)
a%heap_size = a%heap_size-1 a%heap_size = a%heap_size-1_pInt
call heapify(a,1) call heapify(a,1_pInt)
return return
else else
write (*,*) 'PQ_EXTRACT_MAX: error, attempted to pop non-positive PQ' write (*,*) 'PQ_EXTRACT_MAX: error, attempted to pop non-positive PQ'
@ -682,7 +682,7 @@ contains
forall (j=1:tp%n) forall (j=1:tp%n)
tp%ind(j) = j tp%ind(j) = j
end forall end forall
tp%root => build_tree_for_range(tp,1,tp%n, dummy) tp%root => build_tree_for_range(tp,1_pInt,tp%n, dummy)
end subroutine build_tree end subroutine build_tree
recursive function build_tree_for_range(tp,l,u,parent) result (res) recursive function build_tree_for_range(tp,l,u,parent) result (res)
@ -733,7 +733,7 @@ contains
call spread_in_coordinate(tp,i,l,u,res%box(i)) call spread_in_coordinate(tp,i,l,u,res%box(i))
end do end do
res%cut_dim = 0 res%cut_dim = 0
res%cut_val = 0.0 res%cut_val = 0.0_pReal
res%l = l res%l = l
res%u = u res%u = u
res%left =>null() res%left =>null()
@ -749,7 +749,7 @@ contains
! has only a very small difference. This box is not used ! has only a very small difference. This box is not used
! for searching but only for deciding which coordinate to split on. ! for searching but only for deciding which coordinate to split on.
! !
do i=1,dimen do i=1_pInt,dimen
recompute=.true. recompute=.true.
if (associated(parent)) then if (associated(parent)) then
if (i .ne. parent%cut_dim) then if (i .ne. parent%cut_dim) then
@ -771,7 +771,7 @@ contains
if (.false.) then if (.false.) then
! select exact median to have fully balanced tree. ! select exact median to have fully balanced tree.
m = (l+u)/2 m = (l+u)/2_pInt
call select_on_coordinate(tp%the_data,tp%ind,c,m,l,u) call select_on_coordinate(tp%the_data,tp%ind,c,m,l,u)
else else
! !
@ -781,9 +781,9 @@ contains
! !
if (.true.) then if (.true.) then
! actually compute average ! actually compute average
average = sum(tp%the_data(c,tp%ind(l:u))) / real(u-l+1,pReal) average = sum(tp%the_data(c,tp%ind(l:u))) / real(u-l+1_pInt,pReal)
else else
average = (res%box(c)%upper + res%box(c)%lower)/2.0 average = (res%box(c)%upper + res%box(c)%lower)/2.0_pReal
endif endif
res%cut_val = average res%cut_val = average
@ -797,7 +797,7 @@ contains
! res%cut_val = tp%the_data(c,tp%ind(m)) ! res%cut_val = tp%the_data(c,tp%ind(m))
res%left => build_tree_for_range(tp,l,m,res) res%left => build_tree_for_range(tp,l,m,res)
res%right => build_tree_for_range(tp,m+1,u,res) res%right => build_tree_for_range(tp,m+1_pInt,u,res)
if (associated(res%right) .eqv. .false.) then if (associated(res%right) .eqv. .false.) then
res%box = res%left%box res%box = res%left%box
@ -1019,7 +1019,7 @@ contains
type(kdtree2_result), target :: results(:) type(kdtree2_result), target :: results(:)
sr%ballsize = huge(1.0) sr%ballsize = huge(1.0_pReal)
sr%qv => qv sr%qv => qv
sr%nn = nn sr%nn = nn
sr%nfound = 0 sr%nfound = 0
@ -1062,7 +1062,7 @@ contains
allocate (sr%qv(tp%dimen)) allocate (sr%qv(tp%dimen))
sr%qv = tp%the_data(:,idxin) ! copy the vector sr%qv = tp%the_data(:,idxin) ! copy the vector
sr%ballsize = huge(1.0) ! the largest real(pReal) number sr%ballsize = huge(1.0_pReal) ! the largest real(pReal) number
sr%centeridx = idxin sr%centeridx = idxin
sr%correltime = correltime sr%correltime = correltime
@ -1438,7 +1438,7 @@ contains
res = (amin-x)**2; res = (amin-x)**2;
return return
else else
res = 0.0 res = 0.0_pReal
return return
endif endif
endif endif
@ -1461,9 +1461,9 @@ contains
dimen = sr%dimen dimen = sr%dimen
ballsize = sr%ballsize ballsize = sr%ballsize
dis = 0.0 dis = 0.0_pReal
res = .true. res = .true.
do i=1,dimen do i=1_pInt,dimen
l = node%box(i)%lower l = node%box(i)%lower
u = node%box(i)%upper u = node%box(i)%upper
dis = dis + (dis2_from_bnd(sr%qv(i),l,u)) dis = dis + (dis2_from_bnd(sr%qv(i),l,u))
@ -1515,22 +1515,22 @@ contains
mainloop: do i = node%l, node%u mainloop: do i = node%l, node%u
if (rearrange) then if (rearrange) then
sd = 0.0 sd = 0.0_pReal
do k = 1,dimen do k = 1_pInt,dimen
sd = sd + (data(k,i) - qv(k))**2 sd = sd + (data(k,i) - qv(k))**2.0_pReal
if (sd>ballsize) cycle mainloop if (sd>ballsize) cycle mainloop
end do end do
indexofi = ind(i) ! only read it if we have not broken out indexofi = ind(i) ! only read it if we have not broken out
else else
indexofi = ind(i) indexofi = ind(i)
sd = 0.0 sd = 0.0_pReal
do k = 1,dimen do k = 1_pInt,dimen
sd = sd + (data(k,indexofi) - qv(k))**2 sd = sd + (data(k,indexofi) - qv(k))**2.0_pReal
if (sd>ballsize) cycle mainloop if (sd>ballsize) cycle mainloop
end do end do
endif endif
if (centeridx > 0) then ! doing correlation interval? if (centeridx > 0_pInt) then ! doing correlation interval?
if (abs(indexofi-centeridx) < correltime) cycle mainloop if (abs(indexofi-centeridx) < correltime) cycle mainloop
endif endif
@ -1638,26 +1638,26 @@ contains
! which index to the point do we use? ! which index to the point do we use?
if (rearrange) then if (rearrange) then
sd = 0.0 sd = 0.0_pReal
do k = 1,dimen do k = 1_pInt,dimen
sd = sd + (data(k,i) - qv(k))**2 sd = sd + (data(k,i) - qv(k))**2.0_pReal
if (sd>ballsize) cycle mainloop if (sd>ballsize) cycle mainloop
end do end do
indexofi = ind(i) ! only read it if we have not broken out indexofi = ind(i) ! only read it if we have not broken out
else else
indexofi = ind(i) indexofi = ind(i)
sd = 0.0 sd = 0.0_pReal
do k = 1,dimen do k = 1_pInt,dimen
sd = sd + (data(k,indexofi) - qv(k))**2 sd = sd + (data(k,indexofi) - qv(k))**2.0_pReal
if (sd>ballsize) cycle mainloop if (sd>ballsize) cycle mainloop
end do end do
endif endif
if (centeridx > 0) then ! doing correlation interval? if (centeridx > 0_pInt) then ! doing correlation interval?
if (abs(indexofi-centeridx)<correltime) cycle mainloop if (abs(indexofi-centeridx)<correltime) cycle mainloop
endif endif
nfound = nfound+1 nfound = nfound+1_pInt
if (nfound .gt. sr%nalloc) then if (nfound .gt. sr%nalloc) then
! oh nuts, we have to add another one to the tree but ! oh nuts, we have to add another one to the tree but
! there isn't enough room. ! there isn't enough room.
@ -1687,22 +1687,22 @@ contains
real(pReal), allocatable :: all_distances(:) real(pReal), allocatable :: all_distances(:)
! .. ! ..
allocate (all_distances(tp%n)) allocate (all_distances(tp%n))
do i = 1, tp%n do i = 1_pInt, tp%n
all_distances(i) = square_distance(tp%dimen,qv,tp%the_data(:,i)) all_distances(i) = square_distance(tp%dimen,qv,tp%the_data(:,i))
end do end do
! now find 'n' smallest distances ! now find 'n' smallest distances
do i = 1, nn do i = 1_pInt, nn
results(i)%dis = huge(1.0) results(i)%dis = huge(1.0_pReal)
results(i)%idx = -1 results(i)%idx = -1_pInt
end do end do
do i = 1, tp%n do i = 1_pInt, tp%n
if (all_distances(i)<results(nn)%dis) then if (all_distances(i)<results(nn)%dis) then
! insert it somewhere on the list ! insert it somewhere on the list
do j = 1, nn do j = 1, nn
if (all_distances(i)<results(j)%dis) exit if (all_distances(i)<results(j)%dis) exit
end do end do
! now we know 'j' ! now we know 'j'
do k = nn - 1, j, -1 do k = nn - 1_pInt, j, -1_pInt
results(k+1) = results(k) results(k+1) = results(k)
end do end do
results(j)%dis = all_distances(i) results(j)%dis = all_distances(i)