Binary Search

Binary search is a search algorithm that finds the position of a target value within a sorted array. The following binary search routine is taken from NASTRAN, a venerable finite element analysis program from NASA:

      SUBROUTINE BISLOC (*,ID,ARR,LEN,KN,JLOC)
C-----
C     BINARY SEARCH - LOCATE KEY WORD 'ID' IN ARRAY 'ARR', 1ST ENTRY
C     IF FOUND, 'JLOC' IS THE MATCHED POSITION IN 'ARR'
C     IF NOT FOUND, NON-STANDARD RETURN
C                                                     I.E.
C     ID  = KEY WORD TO MATCH IN ARR.      MATCH AGAINST 1ST COL OF ARR
C     ARR = ARRAY TO SEARCH.                          ARR(ROW,COL)
C     LEN = LENGTH OF EACH ENTRY IN ARRAY.            LEN=ROW
C     KN  = NUMBER OF ENTRIES IN THE ARR.             KN =COL
C     JLOC= POINTER RETURNED - FIRST WORD OF ENTRY.   MATCHED ROW
C-----
C
      INTEGER  ARR(1)
      DATA     ISWTCH / 16 /
C
      JJ = LEN - 1
      IF (KN .LT. ISWTCH) GO TO 120
      KLO = 1
      KHI = KN
   10 K   = (KLO+KHI+1)/2
   20 J   = K*LEN - JJ
      IF (ID-ARR(J)) 30,90,40
   30 KHI = K
      GO TO 50
   40 KLO = K
   50 IF (KHI-KLO -1) 100,60,10
   60 IF (K .EQ. KLO) GO TO 70
      K   = KLO
      GO TO 80
   70 K   = KHI
   80 KLO = KHI
      GO TO 20
   90 JLOC = J
      RETURN
  100 JLOC = KHI*LEN - JJ
      J    = KN *LEN - JJ
      IF (ID .GT.ARR(J)) JLOC = JLOC + LEN
  110 RETURN 1
C
C     SEQUENTIAL SEARCH MORE EFFICIENT
C
  120 KHI = KN*LEN - JJ
      DO 130 J = 1,KHI,LEN
      IF (ARR(J)-ID) 130,90,140
  130 CONTINUE
      JLOC = KHI + LEN
      GO TO 110
  140 JLOC = J
      GO TO 110
      END

As you can see, it is a scene of FORTRAN 66 horrors! However, it can be refactored into a nice modern Fortran routine like so:

pure function bisloc(id,arr) result(jloc)

!! binary search of a sorted array.

implicit none

integer,intent(in) :: id
   !! key word to match in `arr`
integer,dimension(:),intent(in) :: arr
   !! array to search (it is
   !! assumed to be sorted)
integer :: jloc
   !! the first matched index in 'arr'
   !! (if not found, 0 is returned)

integer :: j,k,khi,klo,n
integer,parameter :: iswtch = 16

n = size(arr)
jloc = 0

if ( n<iswtch ) then

    ! sequential search more efficient
    
    do j = 1 , n
        if ( arr(j)==id ) then
            jloc = j
            return
        elseif (arr(j)>id) then
            return ! error
        endif
    enddo
    return ! error

else

    klo = 1
    khi = n
    k = (klo+khi+1)/2
    do
        j = k
        if ( id<arr(j) ) then
            khi = k
        elseif ( id==arr(j) ) then
            jloc = j
            return
        else
            klo = k
        endif
        if ( khi-klo<1 ) then
            return ! error
        elseif ( khi-klo==1 ) then
            if ( k==klo ) then
                k = khi
            else
                k = klo
            endif
            klo = khi
        else
            k = (klo+khi+1)/2
        endif
    end do

endif

end function bisloc

I did change a couple of things:

  • The original routine was using an obsolete alternate return (note the * argument) for an error condition. The new one just returns 0 if there is an error.
  • The original routine was passing in the array as a vector and treating it like a matrix (it’s either searching a row or column of this matrix, I admit I had a hard time telling). The new version assumes the input is a rank-1 array. You can pass in whatever matrix slices you want (say a column: a(:,1) or a row: a(1,:)).
  • I also made it a pure function, because that’s how I roll.

Note that NASTRAN was released under the NASA Open Source Agreement version 1.3, a somewhat unusual open source license. It isn’t considered a “free” license by the FSF since it includes a statement requiring changes to be your “original creation”.

Unfortunately, the trend at NASA (and elsewhere) is to declare old Fortran code like this to be totally unmaintainable. So, instead of refactoring them into modern Fortran standards, they are thrown away and rewritten from scratch using a more “modern” language like C++ (whose roots go back to 1972!).

See also

Posted in Programming Tagged with: , , , ,

Leave a Reply

Your email address will not be published. Required fields are marked *

*