C   INCLUDE 'calpgm.inc'
C
C   Herbert M. Pickett, 24 Feb 1989
C   Adapted from subset of BLAS routines in LINPAK 
C
      FUNCTION IDAMAX(N,SX,INCX)
C
C     FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE.
C
      REAL*8  SX(0:*),SMAX
      INTEGER IDAMAX,N,INCX,I,M
      INTEGER*4 IX
C     EMA SX
C
      IDAMAX = 0
      M=N-1
      IF(M.LE.0) THEN
         IF(M.LT.0) IDAMAX=-1
      ELSE IF(INCX.NE.1) THEN
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
         SMAX = ABS(SX(0))
         IX = INCX
         DO 10 I = 1,M
             IF(ABS(SX(IX)).GT.SMAX) THEN
               IDAMAX = I
               SMAX = ABS(SX(IX))
             ENDIF   
             IX = IX + INCX
   10        CONTINUE
      ELSE
C
C        CODE FOR INCREMENT EQUAL TO 1
C
         SMAX = ABS(SX(0))
         DO 30 I = 1,M
            IF(ABS(SX(I)).GT.SMAX) THEN
              IDAMAX = I
              SMAX = ABS(SX(I))
            ENDIF  
   30       CONTINUE
      ENDIF
      IDAMAX=IDAMAX+1
      RETURN
      END
      FUNCTION DASUM(N,SX,INCX)
C
C     TAKES THE SUM OF THE ABSOLUTE VALUES.
C
      REAL*8  DASUM,SX(0:*)
      INTEGER I,INCX,N,M
      INTEGER*4 IX
C     EMA SX
C
      DASUM = 0
      M=N-1
      IF(M.GE.0) THEN
        IF(INCX.NE.1) THEN
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
          IX=0
          DO 10 I = 1,N
            DASUM = DASUM + ABS(SX(IX))
            IX = IX + INCX
   10     CONTINUE
        ELSE
C
C        CODE FOR INCREMENT EQUAL TO 1
C
C       
          DO 30 I = 0,M
            DASUM = DASUM + ABS(SX(I))
   30     CONTINUE
        ENDIF
      ENDIF  
      RETURN
      END
      SUBROUTINE DAXPY(N,SA,SX,INCX,SY,INCY)
C
C     CONSTANT TIMES A VECTOR PLUS A VECTOR.
C
      REAL*8  SX(0:*),SY(0:*),SA
      INTEGER I,INCX,INCY,N,M
      INTEGER*4 IX,IY
C     EMA SX,SY
C
      M=N-1
      IF(M.LT.0)RETURN
      IF(SA.EQ.0) RETURN
      IF(INCX.NE.1)GO TO 20
      IF(INCY.NE.1)GO TO 20
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
      DO 30 I = 0,M
        SY(I) = SY(I) + SA*SX(I)
   30 CONTINUE
      RETURN
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
   20 IX = 0
      IF(INCX.LT.0) IX = -M*INCX
      IY = 0
      IF(INCY.LT.0) IY = -M*INCY
      DO 10 I = 1,N
        SY(IY) = SY(IY) + SA*SX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
      END
      SUBROUTINE  DCOPY(N,SX,INCX,SY,INCY)
C
C     COPIES A VECTOR, X, TO A VECTOR, Y.
C
      REAL*8  SX(0:*),SY(0:*),SA
      INTEGER I,INCX,INCY,N,M
      INTEGER*4 IX,IY
C     EMA SX,SY
C
      M=N-1
      IF(M.LT.0)RETURN
      IF(INCX.EQ.0) THEN
        SA=SX(0)
        IF(INCY.NE.1) THEN
C
C        CODE FOR FILL AND INCREMENT NOT EQUAL TO 1
C
          IY = 0
          IF(INCY.LT.0) IY = -M*INCY
          DO 10 I = 1,N
            SY(IY) = SA
            IY = IY + INCY
   10       CONTINUE
        ELSE
C
C        CODE FOR FILL AND INCREMENT EQUAL TO 1
C
C
          DO 20 I = 0,M
            SY(I) = SA
   20       CONTINUE
        ENDIF
      ELSE IF(INCX.NE.1.OR.INCY.NE.1) THEN
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
        IX = 0
        IF(INCX.LT.0)IX = -M*INCX
        IY = 0
        IF(INCY.LT.0)IY = -M*INCY
        DO 30 I = 1,N
          SY(IY) = SX(IX)
          IX = IX + INCX
          IY = IY + INCY
   30     CONTINUE
      ELSE
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
        DO 40 I = 0,M
          SY(I) = SX(I)
   40     CONTINUE
      ENDIF
      RETURN
      END
      FUNCTION DDOT(N,SX,INCX,SY,INCY)
C
C     FORMS THE DOT PRODUCT OF TWO VECTORS.
C
      REAL*8  DDOT,SX(0:*),SY(0:*)
      INTEGER I,INCX,INCY,N,M
      INTEGER*4 IX,IY
C     EMA SX,SY
C
      DDOT = 0
      M=N-1
      IF(M.LT.0)RETURN
      IF(INCX.NE.1)GO TO 20
      IF(INCY.NE.1)GO TO 20
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
      DO 30 I = 0,M
        DDOT = DDOT + SX(I)*SY(I)
   30 CONTINUE
      RETURN
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
   20 IX = 0
      IF(INCX.LT.0)IX = -M*INCX
      IY = 0
      IF(INCY.LT.0)IY = -M*INCY
      DO 10 I = 1,N
        DDOT = DDOT + SX(IX)*SY(IY)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
      END
      SUBROUTINE  DSCAL(N,SA,SX,INCX)
C
C     SCALES A VECTOR BY A CONSTANT.
C
      REAL*8  SA,SX(0:*)
      INTEGER I,INCX,N,M
      INTEGER*4 IX
C     EMA SX
C     
      M=N-1
      IF(M.LT.0)  RETURN
      IF(SA.EQ.1) RETURN
      IF(INCX.EQ.1) GO TO 20
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
      IX=0
      DO 10 I = 1,N
        SX(IX) = SA*SX(IX)
        IX = IX + INCX
   10 CONTINUE
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
   20 DO 30 I = 0,M
        SX(I) = SA*SX(I)
   30 CONTINUE
      RETURN
      END
      SUBROUTINE  DSWAP (N,SX,INCX,SY,INCY)
C
C     INTERCHANGES TWO VECTORS.
C
      REAL*8  SX(0:*),SY(0:*),STEMP
      INTEGER I,INCX,INCY,N,M
      INTEGER*4 IX,IY
C     EMA SX,SY
C
      M=N-1
      IF(M.LT.0) RETURN
      IF(INCX.NE.1) GO TO 20
C
C       CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
      DO 30 I = 0,M
        STEMP = SX(I)
        SX(I) = SY(I)
        SY(I) = STEMP
   30 CONTINUE
      RETURN
C
C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
C         TO 1
C
   20 IX = 0
      IF(INCX.LT.0)IX = -M*INCX
      IY = 0
      IF(INCY.LT.0)IY = -M*INCY
      DO 10 I = 1,N
        STEMP = SX(IX)
        SX(IX) = SY(IY)
        SY(IY) = STEMP
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
      END
      SUBROUTINE  DROT (N,SX,INCX,SY,INCY,C,S)
C
C     APPLIES A PLANE ROTATION.
C
      REAL*8  SX(0:*),SY(0:*),C,S,STEMP
      INTEGER INCX,INCY,N,I,M
      INTEGER*4 IX,IY
C     EMA SX,SY
C
      M=N-1
      IF(M.LT.0)RETURN
      IF(INCX.NE.1) GO TO 20
      IF(INCY.NE.1) GO TO 20
C
C       CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
      DO 30 I = 0,M
        STEMP = SX(I)*C + SY(I)*S
        SY(I) = SY(I)*C - SX(I)*S
        SX(I) = STEMP
   30 CONTINUE
C
      RETURN
C
C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
C         TO 1
C
   20 IX = 0
      IF(INCX.LT.0)IX = -M*INCX
      IY = 0
      IF(INCY.LT.0)IY = -M*INCY
      DO 10 I = 1,N
        STEMP  = SX(IX)*C + SY(IY)*S
        SY(IY) = SY(IY)*C - SX(IX)*S
        SX(IX) = STEMP
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
      END