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

- NASTRAN-95 at NASA’s GitHub.

## Leave a Reply