**********************************************************************
C
C     Copyright (C) 1992  Roland W. Freund and Noel M. Nachtigal
C     All rights reserved.
C
C     This code is part of a copyrighted package.  For details, see the
C     file `cpyrit.doc' in the top-level directory.
C
C     *****************************************************************
C     ANY USE OF  THIS CODE CONSTITUTES ACCEPTANCE OF  THE TERMS OF THE
C                             COPYRIGHT NOTICE
C     *****************************************************************
C
C**********************************************************************
C
      SUBROUTINE ZAXPBY (N,ZZ,ZA,ZX,ZB,ZY)
C
C     Purpose:
C     This subroutine computes ZZ = ZA * ZX + ZB * ZY.  Several special
C     cases are handled separately:
C        ZA =  0.0, ZB =  0.0 => ZZ = 0.0
C        ZA =  0.0, ZB =  1.0 => ZZ = ZY  (this is COPY)
C        ZA =  0.0, ZB = -1.0 => ZZ = -ZY
C        ZA =  0.0, ZB =   ZB => ZZ = ZB * ZY  (this is SCAL)
C        ZA =  1.0, ZB =  0.0 => ZZ = ZX  (this is COPY)
C        ZA =  1.0, ZB =  1.0 => ZZ = ZX + ZY
C        ZA =  1.0, ZB = -1.0 => ZZ = ZX - ZY
C        ZA =  1.0, ZB =   ZB => ZZ = ZX + ZB * ZY (this is AXPY)
C        ZA = -1.0, ZB =  0.0 => ZZ = -ZX
C        ZA = -1.0, ZB =  1.0 => ZZ = -ZX + ZY
C        ZA = -1.0, ZB = -1.0 => ZZ = -ZX - ZY
C        ZA = -1.0, ZB =   ZB => ZZ = -ZX + ZB * ZY
C        ZA =   ZA, ZB =  0.0 => ZZ = ZA * ZX  (this is SCAL)
C        ZA =   ZA, ZB =  1.0 => ZZ = ZA * ZX + ZY  (this is AXPY)
C        ZA =   ZA, ZB = -1.0 => ZZ = ZA * ZX - ZY
C        ZA =   ZA, ZB =   ZB => ZZ = ZA * ZX + ZB * ZY
C     ZZ may be the same as ZX or ZY.
C
C     Parameters:
C     N  = the dimension of the vectors (input).
C     ZZ = the vector result (output).
C     ZA = scalar multiplier for ZX (input).
C     ZX = one of the vectors (input).
C     ZB = scalar multiplier for ZY (input).
C     ZY = the other vector (input).
C
C     Noel M. Nachtigal
C     March 23, 1993
C
C**********************************************************************
C
      INTRINSIC DIMAG, DREAL
C
      INTEGER N
      DOUBLE COMPLEX ZA, ZB, ZX(N), ZY(N), ZZ(N)
C
C     Local variables.
C
      INTEGER I
      DOUBLE PRECISION DAI, DAR, DBI, DBR
C
      IF (N.LE.0) RETURN
C
      DAI = DIMAG(ZA)
      DAR = DREAL(ZA)
      DBI = DIMAG(ZB)
      DBR = DREAL(ZB)
      IF ((DAR.EQ.0.0D0).AND.(DAI.EQ.0.0D0)) THEN
         IF ((DBR.EQ.0.0D0).AND.(DBI.EQ.0.0D0)) THEN
C           ZA = 0.0, ZB = 0.0 => ZZ = 0.0.
            DO 10 I = 1, N
               ZZ(I) = (0.0D0,0.0D0)
 10         CONTINUE
         ELSE IF ((DBR.EQ.1.0D0).AND.(DBI.EQ.0.0D0)) THEN
C           ZA = 0.0, ZB = 1.0 => ZZ = ZY (this is COPY).
            DO 20 I = 1, N
               ZZ(I) = ZY(I)
 20         CONTINUE
         ELSE IF ((DBR.EQ.-1.0D0).AND.(DBI.EQ.0.0D0)) THEN
C           ZA = 0.0, ZB = -1.0 => ZZ = -ZY.
            DO 30 I = 1, N
               ZZ(I) = -ZY(I)
 30         CONTINUE
         ELSE
C           ZA = 0.0, ZB = ZB => ZZ = ZB * ZY (this is SCAL).
            DO 40 I = 1, N
               ZZ(I) = ZB * ZY(I)
 40         CONTINUE
         END IF
      ELSE IF ((DAR.EQ.1.0D0).AND.(DAI.EQ.0.0D0)) THEN
         IF ((DBR.EQ.0.0D0).AND.(DBI.EQ.0.0D0)) THEN
C           ZA = 1.0, ZB = 0.0 => ZZ = ZX (this is COPY).
            DO 50 I = 1, N
               ZZ(I) = ZX(I)
 50         CONTINUE
         ELSE IF ((DBR.EQ.1.0D0).AND.(DBI.EQ.0.0D0)) THEN
C           ZA = 1.0, ZB = 1.0 => ZZ = ZX + ZY.
            DO 60 I = 1, N
               ZZ(I) = ZX(I) + ZY(I)
 60         CONTINUE
         ELSE IF ((DBR.EQ.-1.0D0).AND.(DBI.EQ.0.0D0)) THEN
C           ZA = 1.0, ZB = -1.0 => ZZ = ZX - ZY.
            DO 70 I = 1, N
               ZZ(I) = ZX(I) - ZY(I)
 70         CONTINUE
         ELSE
C           ZA = 1.0, ZB = ZB => ZZ = ZX + ZB * ZY (this is AXPY).
            DO 80 I = 1, N
               ZZ(I) = ZX(I) + ZB * ZY(I)
 80         CONTINUE
         END IF
      ELSE IF ((DAR.EQ.-1.0D0).AND.(DAI.EQ.0.0D0)) THEN
         IF ((DBR.EQ.0.0D0).AND.(DBI.EQ.0.0D0)) THEN
C           ZA = -1.0, ZB = 0.0 => ZZ = -ZX
            DO 90 I = 1, N
               ZZ(I) = -ZX(I)
 90         CONTINUE
         ELSE IF ((DBR.EQ.1.0D0).AND.(DBI.EQ.0.0D0)) THEN
C           ZA = -1.0, ZB = 1.0 => ZZ = -ZX + ZY
            DO 100 I = 1, N
               ZZ(I) = -ZX(I) + ZY(I)
 100        CONTINUE
         ELSE IF ((DBR.EQ.-1.0D0).AND.(DBI.EQ.0.0D0)) THEN
C           ZA = -1.0, ZB = -1.0 => ZZ = -ZX - ZY.
            DO 110 I = 1, N
               ZZ(I) = -ZX(I) - ZY(I)
 110        CONTINUE
         ELSE
C           ZA = -1.0, ZB = ZB => ZZ = -ZX + ZB * ZY
            DO 120 I = 1, N
               ZZ(I) = -ZX(I) + ZB * ZY(I)
 120        CONTINUE
         END IF
      ELSE
         IF ((DBR.EQ.0.0D0).AND.(DBI.EQ.0.0D0)) THEN
C           ZA = ZA, ZB = 0.0 => ZZ = ZA * ZX (this is SCAL).
            DO 130 I = 1, N
               ZZ(I) = ZA * ZX(I)
 130        CONTINUE
         ELSE IF ((DBR.EQ.1.0D0).AND.(DBI.EQ.0.0D0)) THEN
C           ZA = ZA, ZB = 1.0 => ZZ = ZA * ZX + ZY (this is AXPY)
            DO 140 I = 1, N
               ZZ(I) = ZA * ZX(I) + ZY(I)
 140        CONTINUE
         ELSE IF ((DBR.EQ.-1.0D0).AND.(DBI.EQ.0.0D0)) THEN
C           ZA = ZA, ZB = -1.0 => ZZ = ZA * ZX - ZY.
            DO 150 I = 1, N
               ZZ(I) = ZA * ZX(I) - ZY(I)
 150        CONTINUE
         ELSE
C           ZA = ZA, ZB = ZB => ZZ = ZA * ZX + ZB * ZY.
            DO 160 I = 1, N
               ZZ(I) = ZA * ZX(I) + ZB * ZY(I)
 160        CONTINUE
         END IF
      END IF
C
      RETURN
      END
C
C**********************************************************************