Implement S/C/D/Z AXPBY

BLAS/SRC/CMakeLists.txt

@@ -29,18 +29,18 @@
 #  Level 1 BLAS
 #---------------------------------------------------------
 
-set(SBLAS1 isamax.f sasum.f saxpy.f scopy.f sdot.f snrm2.f90
+set(SBLAS1 isamax.f sasum.f saxpy.f saxpby.f scopy.f sdot.f snrm2.f90
 	srot.f srotg.f90 sscal.f sswap.f sdsdot.f srotmg.f srotm.f)
 
-set(CBLAS1 scabs1.f scasum.f scnrm2.f90 icamax.f caxpy.f ccopy.f
+set(CBLAS1 scabs1.f scasum.f scnrm2.f90 icamax.f caxpy.f caxpby.f ccopy.f
 	cdotc.f cdotu.f csscal.f crotg.f90 cscal.f cswap.f csrot.f)
 
-set(DBLAS1 idamax.f dasum.f daxpy.f dcopy.f ddot.f dnrm2.f90
+set(DBLAS1 idamax.f dasum.f daxpy.f daxpby.f dcopy.f ddot.f dnrm2.f90
 	drot.f drotg.f90 dscal.f dsdot.f dswap.f drotmg.f drotm.f)
 
 set(DB1AUX sscal.f isamax.f)
 
-set(ZBLAS1 dcabs1.f dzasum.f dznrm2.f90 izamax.f zaxpy.f zcopy.f
+set(ZBLAS1 dcabs1.f dzasum.f dznrm2.f90 izamax.f zaxpy.f zaxpby.f zcopy.f
 	zdotc.f zdotu.f zdscal.f zrotg.f90 zscal.f zswap.f zdrot.f)
 
 set(CB1AUX

BLAS/SRC/Makefile

@@ -69,19 +69,19 @@ all: $(BLASLIB)
 #  Comment out the next 6 definitions if you already have
 #  the Level 1 BLAS.
 #---------------------------------------------------------
-SBLAS1 = isamax.o sasum.o saxpy.o scopy.o sdot.o snrm2.o \
+SBLAS1 = isamax.o sasum.o saxpy.o saxpby.o scopy.o sdot.o snrm2.o \
 	srot.o srotg.o sscal.o sswap.o sdsdot.o srotmg.o srotm.o
 $(SBLAS1): $(FRC)
 
-CBLAS1 = scabs1.o scasum.o scnrm2.o icamax.o caxpy.o ccopy.o \
+CBLAS1 = scabs1.o scasum.o scnrm2.o icamax.o caxpy.o caxpby.o ccopy.o \
 	cdotc.o cdotu.o csscal.o crotg.o cscal.o cswap.o csrot.o
 $(CBLAS1): $(FRC)
 
-DBLAS1 = idamax.o dasum.o daxpy.o dcopy.o ddot.o dnrm2.o \
+DBLAS1 = idamax.o dasum.o daxpy.o daxpby.o dcopy.o ddot.o dnrm2.o \
 	drot.o drotg.o dscal.o dsdot.o dswap.o drotmg.o drotm.o
 $(DBLAS1): $(FRC)
 
-ZBLAS1 = dcabs1.o dzasum.o dznrm2.o izamax.o zaxpy.o zcopy.o \
+ZBLAS1 = dcabs1.o dzasum.o dznrm2.o izamax.o zaxpy.o zaxpby.o zcopy.o \
 	zdotc.o zdotu.o zdscal.o zrotg.o zscal.o zswap.o zdrot.o
 $(ZBLAS1): $(FRC)
 

BLAS/SRC/caxpby.f

@@ -0,0 +1,145 @@
+*> \brief \b CAXPBY
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CAXPBY(N,CA,CX,INCX,CB,CY,INCY)
+*
+*       .. Scalar Arguments ..
+*       COMPLEX CA,CB
+*       INTEGER INCX,INCY,N
+*       ..
+*       .. Array Arguments ..
+*       COMPLEX CX(*),CY(*)
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*>    CAXPBY constant times a vector plus constanttimes a vector.
+*>
+*>    Y = ALPHA * X + BETA * Y
+*>
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          number of elements in input vector(s)
+*> \endverbatim
+*>
+*> \param[in] CA
+*> \verbatim
+*>          CA is COMPLEX
+*>          On entry, CA specifies the scalar alpha.
+*> \endverbatim
+*>
+*> \param[in] CX
+*> \verbatim
+*>          CX is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
+*> \endverbatim
+*>
+*> \param[in] INCX
+*> \verbatim
+*>          INCX is INTEGER
+*>          storage spacing between elements of CX
+*> \endverbatim
+*>
+*> \param[in] CB
+*> \verbatim
+*>          CB is COMPLEX
+*>          On entry, CB specifies the scalar beta.
+*> \endverbatim
+*>
+*> \param[in,out] CY
+*> \verbatim
+*>          CY is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
+*> \endverbatim
+*>
+*> \param[in] INCY
+*> \verbatim
+*>          INCY is INTEGER
+*>          storage spacing between elements of CY
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*> \author Martin Koehler, MPI Magdeburg
+*
+*> \ingroup axpby
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>     jack dongarra, linpack, 3/11/78.
+*>     modified 12/3/93, array(1) declarations changed to array(*)
+*>     modfied  8/23/24, implement the axpby case
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE CAXPBY(N,CA,CX,INCX,CB,CY,INCY)
+*
+*  -- Reference BLAS level1 routine --
+*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      COMPLEX CA, CB
+      INTEGER INCX,INCY,N
+*     ..
+*     .. Array Arguments ..
+      COMPLEX CX(*),CY(*)
+*     ..
+*
+*  =====================================================================
+*
+*     .. Local Scalars ..
+      INTEGER I,IX,IY
+*     ..
+      IF (N.LE.0) RETURN
+      IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
+*
+*        code for both increments equal to 1
+*
+         DO I = 1,N
+            CY(I) = CB*CY(I) + CA*CX(I)
+         END DO
+      ELSE
+*
+*        code for unequal increments or equal increments
+*          not equal to 1
+*
+         IX = 1
+         IY = 1
+         IF (INCX.LT.0) IX = (-N+1)*INCX + 1
+         IF (INCY.LT.0) IY = (-N+1)*INCY + 1
+         DO I = 1,N
+            CY(IY) = CB*CY(IY) + CA*CX(IX)
+            IX = IX + INCX
+            IY = IY + INCY
+         END DO
+      END IF
+*
+      RETURN
+*
+*     End of CAXBPY
+*
+      END

BLAS/SRC/daxpby.f

@@ -0,0 +1,162 @@
+*> \brief \b DAXPBY
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DAXPBY(N,DA,DX,INCX,DB,DY,INCY)
+*
+*       .. Scalar Arguments ..
+*       DOUBLE PRECISION DA,DB
+*       INTEGER INCX,INCY,N
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION DX(*),DY(*)
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*>    DAXPBY constant times a vector plus constanttimes a vector.
+*>
+*>    Y = ALPHA * X + BETA * Y
+*>
+*>    uses unrolled loops for increments equal to one.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          number of elements in input vector(s)
+*> \endverbatim
+*>
+*> \param[in] DA
+*> \verbatim
+*>           DA is DOUBLE PRECISION
+*>           On entry, DA specifies the scalar alpha.
+*> \endverbatim
+*>
+*> \param[in] DX
+*> \verbatim
+*>          DX is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
+*> \endverbatim
+*>
+*> \param[in] INCX
+*> \verbatim
+*>          INCX is INTEGER
+*>          storage spacing between elements of DX
+*> \endverbatim
+*>
+*> \param[in] DB
+*> \verbatim
+*>           DB is DOUBLE PRECISION
+*>           On entry, DB specifies the scalar beta.
+*> \endverbatim
+*>
+*> \param[in,out] DY
+*> \verbatim
+*>          DY is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
+*> \endverbatim
+*>
+*> \param[in] INCY
+*> \verbatim
+*>          INCY is INTEGER
+*>          storage spacing between elements of DY
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*> \author Martin Koehler, MPI Magdeburg 
+*
+*> \ingroup axpby
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>     jack dongarra, linpack, 3/11/78.
+*>     modified 12/3/93, array(1) declarations changed to array(*)
+*>     modfied  8/23/24, implement the axpby case
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE DAXPBY(N,DA,DX,INCX,DB,DY,INCY)
+*
+*  -- Reference BLAS level1 routine --
+*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION DA,DB
+      INTEGER INCX,INCY,N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION DX(*),DY(*)
+*     ..
+*
+*  =====================================================================
+*
+*     .. Local Scalars ..
+      INTEGER I,IX,IY,M,MP1
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC MOD
+*     ..
+      IF (N.LE.0) RETURN
+      IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
+*
+*        code for both increments equal to 1
+*
+*
+*        clean-up loop
+*
+         M = MOD(N,4)
+         IF (M.NE.0) THEN
+            DO I = 1,M
+               DY(I) = DB*DY(I) + DA*DX(I)
+            END DO
+         END IF
+         IF (N.LT.4) RETURN
+         MP1 = M + 1
+         DO I = MP1,N,4
+            DY(I)   = DB*DY(I) + DA*DX(I)
+            DY(I+1) = DB*DY(I+1) + DA*DX(I+1)
+            DY(I+2) = DB*DY(I+2) + DA*DX(I+2)
+            DY(I+3) = DB*DY(I+3) + DA*DX(I+3)
+         END DO
+      ELSE
+*
+*        code for unequal increments or equal increments
+*          not equal to 1
+*
+         IX = 1
+         IY = 1
+         IF (INCX.LT.0) IX = (-N+1)*INCX + 1
+         IF (INCY.LT.0) IY = (-N+1)*INCY + 1
+         DO I = 1,N
+          DY(IY) = DB*DY(IY) + DA*DX(IX)
+          IX = IX + INCX
+          IY = IY + INCY
+         END DO
+      END IF
+      RETURN
+*
+*     End of DAXPBY
+*
+      END