C
C $Id: ijpair.F,v 1.2 1998/07/16 16:40:01 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
C     THIS FILE CONTAINS SUBROUTINES IJMAKE, IJFIND, AND BINLOC.  THEY
C     ALLOW FAST LOCATION OF AN ATOM PAIR (IATM,JATM) IN THE PAIRLIST
C     IPAIR.  AFTER THE PAIRLIST HAS BEEN CREATED, IJMAKE SHOULD BE
C     CALLED TO CREATE THREE ARRAYS THAT ALLOW FAST ACCESS TO THE PAIRLIST.
C     AFTER CALLING IJMAKE, SUBROUTINE IJFIND MAY BE USED TO LOCATE THE
C     POSITION OF (IATM,JATM) IN IPAIR.  BINLOC IS A BINARY SEARCH
C     ROUTINE USED BY IJMAKE AND IJFIND.
C
C     name change variables in order to have the same variable names
C     everywhere:: oldname -> newname:
C     .            ipmod1  -> ip1old
C     .            ipack   -> ipold
C     .            ipaddr  -> ijold
C
      SUBROUTINE IJMAKE
C
C
C
C     CREATES THREE PAIRLIST-SIZE ARRAYS (IP1OLD,IPOLD,IJOLD) THAT
C     ALLOW NEAR DIRECT ACCESS OF THE ADDRESS OF ANY ATOM PAIR (IATM,JATM)
C     STORED IN THE PAIRLIST ARRAY IPAIR.
C
C     AFTER CALLING THIS ROUTINE, THEN GIVEN ANY ATOM PAIR (IATM,JATM)
C     THAT IS KNOWN TO BE IN THE PAIRLIST, THE FOLLOWING PROCEDURE
C     WILL LOCATE THE ADDRESS OF THAT PAIR IN THE ARRAY IPAIR:
C
C     LET NPAIRS = IP1(NATOMS+1)-1
C         IJPACK = (IATM*(IATM-3))/2 + JATM + 1
C          IJMOD = MOD(IJPACK,NPAIRS)
C           IMIN = IP1OLD(IJMOD)
C           IMAX = IP1OLD(IJMOD+1)-1
C
C     SEARCH IPOLD(I) OVER THE RANGE I=IMIN,IMAX.
C
C     WHEN IPOLD(I).EQ.IJPACK THEN IPMOD(I) WILL CONTAIN THE POSITION
C     IN IPAIR WHERE (IATM,JATM) IS STORED.
C
C     THE SEARCH FROM IMIN TO IMAX WILL SPAN AT MOST NRATIO ENTRIES, WHERE
C     NRATIO = ((NATOMS*(NATOMS-1)/2)/NPAIRS
C            = RATIO OF POTENTIAL PAIRS TO ACTUAL PAIRS STORED.
C
C     IF A BINARY SEARCH IS USED, THEN THE ENTRY CAN BE LOCATED IN AT
C     MOST LOG2(NRATIO)+1 STEPS.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"
C
      NPAIRS = IP1(NATOMS+1)-1
      NMAX = (NATOMS*(NATOMS-1))/2 - 1
C
C     LOOP OVER THE POSSIBLE VALUES FOR IJMOD.
C
      ISTORE = 1
      DO 500 IJMOD=1,NPAIRS
C
C       IP1OLD(IJMOD) MARKS THE BEGINNING OF THE IJMOD BLOCK IN THE
C       ARRAYS IPOLD AND IJOLD.
C
        IP1OLD(IJMOD) = ISTORE
C
C       IF ALL THE PAIRS HAVE BEEN STORED THEN WE ARE DONE.
C
        IF(ISTORE.GT.NPAIRS) RETURN
C
C       DETERMINE WHICH COMBINATIONS OF (IATM,JATM) WILL YIELD A VALUE
C       OF IJPACK THAT GIVES RISE TO THE CURRENT IJMOD VALUE.  NWRAP
C       KEEPS TRACK OF THE NUMBER OF TIMES THE TARGET IJPACK VALUE
C       WRAPS OR CYCLES THROUGH NPAIRS.
C
        NWRAP = -1
 100    NWRAP = NWRAP + 1
        N = IJMOD - 1 + NWRAP*NPAIRS
        IF(N.GT.NMAX) GO TO 500
C
C       TARGET IJPACK VALUE:
C
        IJPACK = N + 1
C
C       MUST SOLVE THE FOLLOWING  QUADRATIC TO DETERMINE THE IATM VALUE
C       THAT WILL YIELD THE TARGET IJPACK VALUE:
C
C       (IATM*(IATM-3))/2 + 2 = IJPACK
C
        DISCR = 2*N + 0.25D0
        XIATM = DSQRT(DISCR) + 1.50000001D0
        IATM = XIATM
C
C       SKIP IATM IF IT HAS NO PAIRS IN THE ACTUAL PAIRLIST.
C
        ISTART = IP1(IATM)
        ISTOP = IP1(IATM+1)-1
        IF(ISTOP.LT.ISTART) GO TO 100
C
C       DETERMINE JATM BASED ON THE LOCATION OF THE (IATM,1) PAIR
C       IN THE PACKED STORAGE SCHEME:
C
        NII = (IATM*(IATM-3))/2
        JATM = N - NII
C
C       DO A BINARY SEARCH OF THE IATM SECTION OF THE ACTUAL PAIRLIST TO
C       DETERMINE WHETHER THE (IATM,JATM) PAIR HAS BEEN STORED.  IF THE
C       PAIR HAS BEEN STORED, THEN IADDR WILL CONTAIN THE CORRESPONDING
C       ADDRESS IN IPAIR.
C
        CALL BINLOC(IPAIR,ISTART,ISTOP,JATM,IADDR)
        IF(IADDR.NE.0)THEN
          IJOLD(ISTORE) = IADDR
          IPOLD(ISTORE) = IJPACK
          ISTORE = ISTORE + 1
        ENDIF
        GO TO 100
 500  CONTINUE
C
C     OVERFLOW POINTER IN CASE ALL POSSIBLE PAIRS HAVE BEEN STORED:
C
      IP1OLD(NPAIRS+1) = NPAIRS+1
      RETURN
      END
C
C
C
      SUBROUTINE IJFIND(NPAIRS,IATM,JATM,IJADDR)
C
C     LOCATES THE ATOM PAIR (IATM,JATM) IN THE PAIRLIST AFTER THE ARRAYS
C     IP1OLD, IPOLD, AND IJOLD HAVE BEEN CREATED BY SUBROUTINE IJMAKE.
C
C     INPUT:
C
C     NPAIRS = TOTAL NUMBER OF STORED PAIRS.
C
C     IATM,
C     JATM   = ATOM PAIR TO BE LOCATED (IATM > JATM).
C
C     RETURNED:
C
C     IJADDR = POSITION IN IPAIR WHERE (IATM,JATM) IS STORED.
C
C     IF (IATM,JATM) IS NOT STORED, THEN IJADDR IS RETURNED WITH
C     A VALUE OF ZERO.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"
C
C     INDEX FOR PACKED, LOWER TRIAGULAR STORAGE:
C
      IJPACK = (IATM*(IATM-3))/2 + JATM + 1
C
C     REMAINDER OF IJPACK/NPAIRS:
C
      IJMOD = MOD(IJPACK,NPAIRS)
      IF(IJMOD.EQ.0) IJMOD = NPAIRS
C
C     DO A BINARY SEARCH IN THE ARRAY IPOLD TO FIND A MATCHING VALUE
C     OF IJPACK.
C
      IMIN = IP1OLD(IJMOD)
      IMAX = IP1OLD(IJMOD+1) - 1
      CALL BINLOC(IPOLD,IMIN,IMAX,IJPACK,IMATCH)
      IF(IMATCH.NE.0)THEN
        IJADDR = IJOLD(IMATCH)
      ELSE
        IJADDR = 0
      ENDIF
      RETURN
      END
C
C
C
      SUBROUTINE BINLOC(LIST,ISTART,ISTOP,ITARGT,IFOUND)
C
C     USES A BINARY SEARCH TO LOCATE A TARGET ELEMENT IN A SORTED LIST.
C
C     LIST   = SORTED LIST OF INTEGERS TO BE SEARCHED.
C     ISTART = STARTING POINT OF SEARCH IN LIST.
C     ISTOP  = STOPPING POINT OF SEARCH IN LIST.
C     ITARGT = VALUE OF TARGET ELEMENT SOUGHT IN LIST.
C     IFOUND = RETURNED LOCATION OF ITARGET IN LIST.  IF NOT PRESENT, THEN
C              IFOUND IS RETURNED WITH A VALUE OF 0.
C
      DIMENSION LIST(*)
C
      IFOUND = 0
C
C     USE BINARY SEARCH ONLY IF LIST HAS MORE THAN ONE VALUE.
C
      IF(ISTART.EQ.ISTOP)THEN
        IF(LIST(ISTART).EQ.ITARGT) IFOUND = ISTART
      ELSE
        IBEGIN = ISTART
        IEND = ISTOP
 10     IHALF = IBEGIN + (IEND-IBEGIN)/2
        LHALF = LIST(IHALF)
        IF(LHALF.EQ.ITARGT)THEN
          IFOUND = IHALF
          RETURN
        ELSEIF(LHALF.GT.ITARGT)THEN
          IEND = IHALF - 1
        ELSE
          IBEGIN = IHALF + 1
        ENDIF
        IF(IBEGIN.GT.IEND)THEN
          RETURN
        ELSE
          GO TO 10
        ENDIF
      ENDIF
      RETURN
      END