DOUG 0.2

BHeap.f90

Go to the documentation of this file.
00001 ! DOUG - Domain decomposition On Unstructured Grids
00002 ! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
00003 ! Department of Mathematics, University of Bath
00004 !
00005 ! This library is free software; you can redistribute it and/or
00006 ! modify it under the terms of the GNU Lesser General Public
00007 ! License as published by the Free Software Foundation; either
00008 ! version 2.1 of the License, or (at your option) any later version.
00009 !
00010 ! This library is distributed in the hope that it will be useful,
00011 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013 ! Lesser General Public License for more details.
00014 !
00015 ! You should have received a copy of the GNU Lesser General Public
00016 ! License along with this library; if not, write to the Free Software
00017 ! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00018 ! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
00019 ! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
00020 ! mailto:info(at)dougdevel.org)
00021 
00022 module BinaryHeap
00023    
00024     implicit none
00025 
00026     !! Max-first binary heap
00027     type BHeap
00028         private
00029         integer :: size ! The number of elements in the structure
00030         integer, pointer :: by(:) ! Numbers by which to sort
00031         integer, pointer :: inds(:) ! Indices to sort
00032     end type
00033 
00034     private :: BHeap_Heapify
00035 
00036 contains
00037 
00038   function BHeap_new() result(B)
00039     type(BHeap) :: B
00040     B%size = 0
00041     B%by => NULL()
00042     B%inds => NULL()
00043   end function BHeap_new
00044 
00045     !! Allocate the structure (initially empty)
00046     subroutine BHeap_init(B, max)
00047         implicit none
00048         
00049         !! The BHeap to allocate
00050         type(BHeap), intent(inout) :: B
00051         !! How many nodes it will have room for
00052         integer, intent(in) :: max
00053         
00054         if (associated(B%by)) deallocate(B%by)
00055         if (associated(B%inds)) deallocate(B%inds)
00056 
00057         B%size=0
00058         
00059         allocate(B%by(max),B%inds(max))
00060                
00061     end subroutine BHeap_init
00062 
00063     !! Destroy the structure
00064     subroutine BHeap_destroy(B)
00065         implicit none
00066 
00067         !! The BHeap to destroy
00068         type(BHeap), intent(inout) :: B
00069 
00070         B%size=-1
00071         if (associated(B%by)) deallocate(B%by)
00072         if (associated(B%inds)) deallocate(B%inds)
00073         
00074     end subroutine BHeap_destroy
00075 
00076     !! Get the max value from the heap
00077     function BHeap_maxv(B) result (res)
00078         implicit none
00079         !! The BHeap to get the value from
00080         type(BHeap), intent(in) :: B
00081         !! Maximum value currently in the heap
00082         integer :: res
00083 
00084         res=B%by(1)
00085     end function BHeap_maxv
00086     
00087     !! Get the index corresponding to max value from the heap
00088     function BHeap_maxi(B) result (res)
00089         implicit none
00090 
00091         !! The BHeap to get the index from
00092         type(BHeap), intent(in) :: B
00093         !! Index corresponding to the maximum value
00094         integer :: res
00095 
00096         res=B%inds(1)
00097     end function BHeap_maxi
00098         
00099     !! Get the current heap size
00100     function BHeap_size(B) result (res)
00101         implicit none
00102 
00103         !! The BHeap to get the size of
00104         type(BHeap), intent(in) :: B
00105         !! Current heap size
00106         integer :: res
00107 
00108         res=B%size
00109     end function BHeap_size
00110 
00111     !! Check if the heap is full
00112     function BHeap_full(B) result (res)
00113         implicit none
00114 
00115         !! The BHeap to check
00116         type(BHeap), intent(in) :: B
00117         logical :: res
00118 
00119         res=(B%size < ubound(B%by,1))
00120     end function BHeap_full
00121     
00122     !! Insert a node into the heap (if there is room)
00123     subroutine BHeap_insert(B, val, ind)
00124         use RealKind
00125         
00126         implicit none
00127 
00128         !! The BHeap to insert to
00129         type(BHeap), intent(inout) :: B
00130         !! The value of the key to be inserted
00131         integer, intent(in) :: val
00132         !! Index associated with the value to be inserted
00133         integer, intent(in) :: ind
00134 
00135         integer :: cur
00136 
00137         if (B%size<ubound(B%by,1)) then
00138             B%size=B%size+1; cur=B%size
00139         
00140             ! Float the value downward
00141             do while (cur>1)
00142                 if (B%by(cur/2)>=val) exit
00143                 B%by(cur)=B%by(cur/2)
00144                 B%inds(cur)=B%inds(cur/2)
00145                 cur=cur/2
00146             enddo
00147         
00148             ! And put the value where appropriate
00149             B%by(cur)=val; B%inds(cur)=ind
00150         endif
00151     end subroutine BHeap_insert
00152 
00153     !! Float element with index i up
00154     subroutine BHeap_heapify(B, curi)
00155          use RealKind
00156         
00157         implicit none
00158 
00159         !! BHeap to fix up
00160         type(BHeap), intent(inout) :: B
00161         !! Index at which to start realing
00162         integer, intent(in) :: curi
00163         
00164         integer :: val
00165         integer :: ind, cur
00166 
00167         integer :: l, r, bt ! left, right, best
00168 
00169         cur=curi
00170 
00171         val=B%by(cur); ind=B%inds(cur)
00172 
00173         bt=cur ! since for every other iteration it is guaranteed
00174         
00175         do
00176             l=2*cur; r=2*cur+1 
00177             ! Find which is the largest of the three
00178             if (r<=B%size) then
00179                  if ( B%by(bt)<B%by(l) ) bt=l
00180                  if ( B%by(bt)<B%by(r) ) bt=r
00181             elseif (l<B%size) then
00182                 if (B%by(bt)<B%by(l)) bt=l
00183             endif
00184 
00185             ! If current is the largest, we have a heap
00186             if (bt==cur) exit
00187          
00188             ! Otherwise move the best down
00189             B%by(cur)=B%by(bt)
00190             B%inds(cur)=B%inds(bt)
00191             cur=bt
00192             
00193             ! And try to place our value to the new position
00194             B%by(cur)=val
00195         enddo
00196         
00197         ! And use the element
00198         B%inds(cur)=ind
00199     
00200     end subroutine BHeap_heapify
00201 
00202     !! Create heap from two arrays
00203     subroutine BHeap_create(B, byvals, inds, size)
00204         use RealKind
00205         
00206         implicit none
00207         
00208         !! The array of values
00209         integer, intent(in) :: byvals(:)
00210         !! The array of indices
00211         integer, intent(in) :: inds(:)
00212         !! The size of the heap to be created
00213         integer, intent(in) :: size
00214         type(BHeap), intent(inout) :: B
00215 
00216         integer :: i
00217         
00218         ! Init the heap with the elements given
00219         call BHeap_init(B,size)
00220         B%size=size
00221         B%by=byvals; B%inds=inds
00222 
00223         ! And make it into a proper heap
00224         do i=B%size/2,1,-1
00225             call BHeap_heapify(B,i)
00226         enddo
00227     end subroutine BHeap_create
00228 
00229     !! Remove the maximal element from the heap
00230     subroutine BHeap_delmax(B)
00231         implicit none
00232 
00233         !! The BHeap to remove max element from
00234         type(BHeap), intent(inout) :: B
00235 
00236         if (B%size>=1) then
00237             ! Move the last element to first
00238             B%by(1)=B%by(B%size); B%inds(1)=B%inds(B%size)
00239             B%size=B%size-1
00240             ! And reel it upwards
00241             call BHeap_heapify(B,1)
00242         endif
00243     end subroutine BHeap_delmax
00244 
00245     !! Heapsort (inplace) - sort both arrays by values array
00246     subroutine HeapSort(byvals, inds, size)
00247         use RealKind
00248         
00249         implicit none
00250 
00251         !! The values array
00252         integer, intent(inout), target :: byvals(:)
00253         !! The indices array
00254         integer, intent(inout), target :: inds(:)
00255         !! The size of the previous two (sort this much)
00256         integer, intent(in) :: size
00257         
00258         type(BHeap)  :: B
00259         integer :: i, ind
00260         real(kind=xyzk) :: val
00261         
00262         ! Create a proper heap
00263         !call BHeap_create(B,byvals,inds, size)
00264 
00265         ! Use byvals and inds as BHeap internal arrays - avoids copying data
00266         B%size=size
00267         B%by=>byvals
00268         B%inds=>inds
00269         
00270         ! And make it into a proper heap
00271         do i=B%size/2,1,-1
00272             call BHeap_heapify(B,i)
00273         enddo
00274 
00275         ! And use it to sort the array 
00276         do i=size,2,-1
00277             ! Remmember the max
00278             val=B%by(1); ind=B%inds(1)
00279             ! Remove it from the heap, freeing up space in the end
00280             call BHeap_delmax(B)
00281             ! And put that data in the end
00282             byvals(i)=val; inds(i)=ind
00283         enddo
00284     end subroutine HeapSort
00285 
00286 end module BinaryHeap