Fix rank-1 update for self-adjoint packed matrices.
diff --git a/blas/CMakeLists.txt b/blas/CMakeLists.txt
index 3877e12..c35a2fd 100644
--- a/blas/CMakeLists.txt
+++ b/blas/CMakeLists.txt
@@ -18,9 +18,9 @@
set(EigenBlas_SRCS ${EigenBlas_SRCS}
complexdots.f
srotm.f srotmg.f drotm.f drotmg.f
- lsame.f dspmv.f ssbmv.f
- chbmv.f chpr.f sspmv.f
- zhbmv.f zhpr.f chpmv.f dsbmv.f
+ lsame.f dspmv.f ssbmv.f
+ chbmv.f sspmv.f
+ zhbmv.f chpmv.f dsbmv.f
zhpmv.f
dtbmv.f stbmv.f ctbmv.f ztbmv.f
)
diff --git a/blas/PackedSelfadjointProduct.h b/blas/PackedSelfadjointProduct.h
index 1ba67b9..f7c9b93 100644
--- a/blas/PackedSelfadjointProduct.h
+++ b/blas/PackedSelfadjointProduct.h
@@ -14,12 +14,6 @@
/* Optimized matrix += alpha * uv'
* The matrix is in packed form.
- *
- * FIXME I always fail tests for complex self-adjoint matrices.
- *
- * ******* FATAL ERROR - PARAMETER NUMBER 6 WAS CHANGED INCORRECTLY *******
- * ******* xHPR FAILED ON CALL NUMBER:
- * 2: xHPR ('U', 1, 0.0, X, 1, AP)
*/
template<typename Scalar, typename Index, int StorageOrder, int UpLo, bool ConjLhs, bool ConjRhs>
struct selfadjoint_packed_rank1_update;
@@ -27,20 +21,20 @@
template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
struct selfadjoint_packed_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRhs>
{
- static void run(Index size, Scalar* mat, const Scalar* vec, Scalar alpha)
+ typedef typename NumTraits<Scalar>::Real RealScalar;
+ static void run(Index size, Scalar* mat, const Scalar* vec, RealScalar alpha)
{
typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
typedef typename conj_expr_if<ConjLhs,OtherMap>::type ConjRhsType;
conj_if<ConjRhs> cj;
- Index offset = 0;
for (Index i=0; i<size; ++i)
{
- Map<Matrix<Scalar,Dynamic,1> >(mat+offset, UpLo==Lower ? size-i : (i+1))
+ Map<Matrix<Scalar,Dynamic,1> >(mat, UpLo==Lower ? size-i : (i+1))
+= alpha * cj(vec[i]) * ConjRhsType(OtherMap(vec+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1)));
//FIXME This should be handled outside.
- mat[offset+(UpLo==Lower ? 0 : i)] = real(mat[offset+(UpLo==Lower ? 0 : i)]);
- offset += UpLo==Lower ? size-i : (i+1);
+ mat[UpLo==Lower ? 0 : i] = real(mat[UpLo==Lower ? 0 : i]);
+ mat += UpLo==Lower ? size-i : (i+1);
}
}
};
@@ -48,7 +42,8 @@
template<typename Scalar, typename Index, int UpLo, bool ConjLhs, bool ConjRhs>
struct selfadjoint_packed_rank1_update<Scalar,Index,RowMajor,UpLo,ConjLhs,ConjRhs>
{
- static void run(Index size, Scalar* mat, const Scalar* vec, Scalar alpha)
+ typedef typename NumTraits<Scalar>::Real RealScalar;
+ static void run(Index size, Scalar* mat, const Scalar* vec, RealScalar alpha)
{
selfadjoint_packed_rank1_update<Scalar,Index,ColMajor,UpLo==Lower?Upper:Lower,ConjRhs,ConjLhs>::run(size,mat,vec,alpha);
}
diff --git a/blas/chpr.f b/blas/chpr.f
deleted file mode 100644
index 11bd5c6..0000000
--- a/blas/chpr.f
+++ /dev/null
@@ -1,220 +0,0 @@
- SUBROUTINE CHPR(UPLO,N,ALPHA,X,INCX,AP)
-* .. Scalar Arguments ..
- REAL ALPHA
- INTEGER INCX,N
- CHARACTER UPLO
-* ..
-* .. Array Arguments ..
- COMPLEX AP(*),X(*)
-* ..
-*
-* Purpose
-* =======
-*
-* CHPR performs the hermitian rank 1 operation
-*
-* A := alpha*x*conjg( x' ) + A,
-*
-* where alpha is a real scalar, x is an n element vector and A is an
-* n by n hermitian matrix, supplied in packed form.
-*
-* Arguments
-* ==========
-*
-* UPLO - CHARACTER*1.
-* On entry, UPLO specifies whether the upper or lower
-* triangular part of the matrix A is supplied in the packed
-* array AP as follows:
-*
-* UPLO = 'U' or 'u' The upper triangular part of A is
-* supplied in AP.
-*
-* UPLO = 'L' or 'l' The lower triangular part of A is
-* supplied in AP.
-*
-* Unchanged on exit.
-*
-* N - INTEGER.
-* On entry, N specifies the order of the matrix A.
-* N must be at least zero.
-* Unchanged on exit.
-*
-* ALPHA - REAL .
-* On entry, ALPHA specifies the scalar alpha.
-* Unchanged on exit.
-*
-* X - COMPLEX array of dimension at least
-* ( 1 + ( n - 1 )*abs( INCX ) ).
-* Before entry, the incremented array X must contain the n
-* element vector x.
-* Unchanged on exit.
-*
-* INCX - INTEGER.
-* On entry, INCX specifies the increment for the elements of
-* X. INCX must not be zero.
-* Unchanged on exit.
-*
-* AP - COMPLEX array of DIMENSION at least
-* ( ( n*( n + 1 ) )/2 ).
-* Before entry with UPLO = 'U' or 'u', the array AP must
-* contain the upper triangular part of the hermitian matrix
-* packed sequentially, column by column, so that AP( 1 )
-* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
-* and a( 2, 2 ) respectively, and so on. On exit, the array
-* AP is overwritten by the upper triangular part of the
-* updated matrix.
-* Before entry with UPLO = 'L' or 'l', the array AP must
-* contain the lower triangular part of the hermitian matrix
-* packed sequentially, column by column, so that AP( 1 )
-* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
-* and a( 3, 1 ) respectively, and so on. On exit, the array
-* AP is overwritten by the lower triangular part of the
-* updated matrix.
-* Note that the imaginary parts of the diagonal elements need
-* not be set, they are assumed to be zero, and on exit they
-* are set to zero.
-*
-* Further Details
-* ===============
-*
-* Level 2 Blas routine.
-*
-* -- Written on 22-October-1986.
-* Jack Dongarra, Argonne National Lab.
-* Jeremy Du Croz, Nag Central Office.
-* Sven Hammarling, Nag Central Office.
-* Richard Hanson, Sandia National Labs.
-*
-* =====================================================================
-*
-* .. Parameters ..
- COMPLEX ZERO
- PARAMETER (ZERO= (0.0E+0,0.0E+0))
-* ..
-* .. Local Scalars ..
- COMPLEX TEMP
- INTEGER I,INFO,IX,J,JX,K,KK,KX
-* ..
-* .. External Functions ..
- LOGICAL LSAME
- EXTERNAL LSAME
-* ..
-* .. External Subroutines ..
- EXTERNAL XERBLA
-* ..
-* .. Intrinsic Functions ..
- INTRINSIC CONJG,REAL
-* ..
-*
-* Test the input parameters.
-*
- INFO = 0
- IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
- INFO = 1
- ELSE IF (N.LT.0) THEN
- INFO = 2
- ELSE IF (INCX.EQ.0) THEN
- INFO = 5
- END IF
- IF (INFO.NE.0) THEN
- CALL XERBLA('CHPR ',INFO)
- RETURN
- END IF
-*
-* Quick return if possible.
-*
- IF ((N.EQ.0) .OR. (ALPHA.EQ.REAL(ZERO))) RETURN
-*
-* Set the start point in X if the increment is not unity.
-*
- IF (INCX.LE.0) THEN
- KX = 1 - (N-1)*INCX
- ELSE IF (INCX.NE.1) THEN
- KX = 1
- END IF
-*
-* Start the operations. In this version the elements of the array AP
-* are accessed sequentially with one pass through AP.
-*
- KK = 1
- IF (LSAME(UPLO,'U')) THEN
-*
-* Form A when upper triangle is stored in AP.
-*
- IF (INCX.EQ.1) THEN
- DO 20 J = 1,N
- IF (X(J).NE.ZERO) THEN
- TEMP = ALPHA*CONJG(X(J))
- K = KK
- DO 10 I = 1,J - 1
- AP(K) = AP(K) + X(I)*TEMP
- K = K + 1
- 10 CONTINUE
- AP(KK+J-1) = REAL(AP(KK+J-1)) + REAL(X(J)*TEMP)
- ELSE
- AP(KK+J-1) = REAL(AP(KK+J-1))
- END IF
- KK = KK + J
- 20 CONTINUE
- ELSE
- JX = KX
- DO 40 J = 1,N
- IF (X(JX).NE.ZERO) THEN
- TEMP = ALPHA*CONJG(X(JX))
- IX = KX
- DO 30 K = KK,KK + J - 2
- AP(K) = AP(K) + X(IX)*TEMP
- IX = IX + INCX
- 30 CONTINUE
- AP(KK+J-1) = REAL(AP(KK+J-1)) + REAL(X(JX)*TEMP)
- ELSE
- AP(KK+J-1) = REAL(AP(KK+J-1))
- END IF
- JX = JX + INCX
- KK = KK + J
- 40 CONTINUE
- END IF
- ELSE
-*
-* Form A when lower triangle is stored in AP.
-*
- IF (INCX.EQ.1) THEN
- DO 60 J = 1,N
- IF (X(J).NE.ZERO) THEN
- TEMP = ALPHA*CONJG(X(J))
- AP(KK) = REAL(AP(KK)) + REAL(TEMP*X(J))
- K = KK + 1
- DO 50 I = J + 1,N
- AP(K) = AP(K) + X(I)*TEMP
- K = K + 1
- 50 CONTINUE
- ELSE
- AP(KK) = REAL(AP(KK))
- END IF
- KK = KK + N - J + 1
- 60 CONTINUE
- ELSE
- JX = KX
- DO 80 J = 1,N
- IF (X(JX).NE.ZERO) THEN
- TEMP = ALPHA*CONJG(X(JX))
- AP(KK) = REAL(AP(KK)) + REAL(TEMP*X(JX))
- IX = JX
- DO 70 K = KK + 1,KK + N - J
- IX = IX + INCX
- AP(K) = AP(K) + X(IX)*TEMP
- 70 CONTINUE
- ELSE
- AP(KK) = REAL(AP(KK))
- END IF
- JX = JX + INCX
- KK = KK + N - J + 1
- 80 CONTINUE
- END IF
- END IF
-*
- RETURN
-*
-* End of CHPR .
-*
- END
diff --git a/blas/dspr.f b/blas/dspr.f
deleted file mode 100644
index 538e4f7..0000000
--- a/blas/dspr.f
+++ /dev/null
@@ -1,202 +0,0 @@
- SUBROUTINE DSPR(UPLO,N,ALPHA,X,INCX,AP)
-* .. Scalar Arguments ..
- DOUBLE PRECISION ALPHA
- INTEGER INCX,N
- CHARACTER UPLO
-* ..
-* .. Array Arguments ..
- DOUBLE PRECISION AP(*),X(*)
-* ..
-*
-* Purpose
-* =======
-*
-* DSPR performs the symmetric rank 1 operation
-*
-* A := alpha*x*x' + A,
-*
-* where alpha is a real scalar, x is an n element vector and A is an
-* n by n symmetric matrix, supplied in packed form.
-*
-* Arguments
-* ==========
-*
-* UPLO - CHARACTER*1.
-* On entry, UPLO specifies whether the upper or lower
-* triangular part of the matrix A is supplied in the packed
-* array AP as follows:
-*
-* UPLO = 'U' or 'u' The upper triangular part of A is
-* supplied in AP.
-*
-* UPLO = 'L' or 'l' The lower triangular part of A is
-* supplied in AP.
-*
-* Unchanged on exit.
-*
-* N - INTEGER.
-* On entry, N specifies the order of the matrix A.
-* N must be at least zero.
-* Unchanged on exit.
-*
-* ALPHA - DOUBLE PRECISION.
-* On entry, ALPHA specifies the scalar alpha.
-* Unchanged on exit.
-*
-* X - DOUBLE PRECISION array of dimension at least
-* ( 1 + ( n - 1 )*abs( INCX ) ).
-* Before entry, the incremented array X must contain the n
-* element vector x.
-* Unchanged on exit.
-*
-* INCX - INTEGER.
-* On entry, INCX specifies the increment for the elements of
-* X. INCX must not be zero.
-* Unchanged on exit.
-*
-* AP - DOUBLE PRECISION array of DIMENSION at least
-* ( ( n*( n + 1 ) )/2 ).
-* Before entry with UPLO = 'U' or 'u', the array AP must
-* contain the upper triangular part of the symmetric matrix
-* packed sequentially, column by column, so that AP( 1 )
-* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
-* and a( 2, 2 ) respectively, and so on. On exit, the array
-* AP is overwritten by the upper triangular part of the
-* updated matrix.
-* Before entry with UPLO = 'L' or 'l', the array AP must
-* contain the lower triangular part of the symmetric matrix
-* packed sequentially, column by column, so that AP( 1 )
-* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
-* and a( 3, 1 ) respectively, and so on. On exit, the array
-* AP is overwritten by the lower triangular part of the
-* updated matrix.
-*
-* Further Details
-* ===============
-*
-* Level 2 Blas routine.
-*
-* -- Written on 22-October-1986.
-* Jack Dongarra, Argonne National Lab.
-* Jeremy Du Croz, Nag Central Office.
-* Sven Hammarling, Nag Central Office.
-* Richard Hanson, Sandia National Labs.
-*
-* =====================================================================
-*
-* .. Parameters ..
- DOUBLE PRECISION ZERO
- PARAMETER (ZERO=0.0D+0)
-* ..
-* .. Local Scalars ..
- DOUBLE PRECISION TEMP
- INTEGER I,INFO,IX,J,JX,K,KK,KX
-* ..
-* .. External Functions ..
- LOGICAL LSAME
- EXTERNAL LSAME
-* ..
-* .. External Subroutines ..
- EXTERNAL XERBLA
-* ..
-*
-* Test the input parameters.
-*
- INFO = 0
- IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
- INFO = 1
- ELSE IF (N.LT.0) THEN
- INFO = 2
- ELSE IF (INCX.EQ.0) THEN
- INFO = 5
- END IF
- IF (INFO.NE.0) THEN
- CALL XERBLA('DSPR ',INFO)
- RETURN
- END IF
-*
-* Quick return if possible.
-*
- IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
-*
-* Set the start point in X if the increment is not unity.
-*
- IF (INCX.LE.0) THEN
- KX = 1 - (N-1)*INCX
- ELSE IF (INCX.NE.1) THEN
- KX = 1
- END IF
-*
-* Start the operations. In this version the elements of the array AP
-* are accessed sequentially with one pass through AP.
-*
- KK = 1
- IF (LSAME(UPLO,'U')) THEN
-*
-* Form A when upper triangle is stored in AP.
-*
- IF (INCX.EQ.1) THEN
- DO 20 J = 1,N
- IF (X(J).NE.ZERO) THEN
- TEMP = ALPHA*X(J)
- K = KK
- DO 10 I = 1,J
- AP(K) = AP(K) + X(I)*TEMP
- K = K + 1
- 10 CONTINUE
- END IF
- KK = KK + J
- 20 CONTINUE
- ELSE
- JX = KX
- DO 40 J = 1,N
- IF (X(JX).NE.ZERO) THEN
- TEMP = ALPHA*X(JX)
- IX = KX
- DO 30 K = KK,KK + J - 1
- AP(K) = AP(K) + X(IX)*TEMP
- IX = IX + INCX
- 30 CONTINUE
- END IF
- JX = JX + INCX
- KK = KK + J
- 40 CONTINUE
- END IF
- ELSE
-*
-* Form A when lower triangle is stored in AP.
-*
- IF (INCX.EQ.1) THEN
- DO 60 J = 1,N
- IF (X(J).NE.ZERO) THEN
- TEMP = ALPHA*X(J)
- K = KK
- DO 50 I = J,N
- AP(K) = AP(K) + X(I)*TEMP
- K = K + 1
- 50 CONTINUE
- END IF
- KK = KK + N - J + 1
- 60 CONTINUE
- ELSE
- JX = KX
- DO 80 J = 1,N
- IF (X(JX).NE.ZERO) THEN
- TEMP = ALPHA*X(JX)
- IX = JX
- DO 70 K = KK,KK + N - J
- AP(K) = AP(K) + X(IX)*TEMP
- IX = IX + INCX
- 70 CONTINUE
- END IF
- JX = JX + INCX
- KK = KK + N - J + 1
- 80 CONTINUE
- END IF
- END IF
-*
- RETURN
-*
-* End of DSPR .
-*
- END
diff --git a/blas/level2_cplx_impl.h b/blas/level2_cplx_impl.h
index 11ee13b..f52d384 100644
--- a/blas/level2_cplx_impl.h
+++ b/blas/level2_cplx_impl.h
@@ -108,10 +108,49 @@
* where alpha is a real scalar, x is an n element vector and A is an
* n by n hermitian matrix, supplied in packed form.
*/
-// int EIGEN_BLAS_FUNC(hpr)(char *uplo, int *n, RealScalar *alpha, RealScalar *x, int *incx, RealScalar *ap)
-// {
-// return 1;
-// }
+int EIGEN_BLAS_FUNC(hpr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pap)
+{
+ typedef void (*functype)(int, Scalar*, const Scalar*, RealScalar);
+ static functype func[2];
+
+ static bool init = false;
+ if(!init)
+ {
+ for(int k=0; k<2; ++k)
+ func[k] = 0;
+
+ func[UP] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run);
+ func[LO] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run);
+
+ init = true;
+ }
+
+ Scalar* x = reinterpret_cast<Scalar*>(px);
+ Scalar* ap = reinterpret_cast<Scalar*>(pap);
+ RealScalar alpha = *palpha;
+
+ int info = 0;
+ if(UPLO(*uplo)==INVALID) info = 1;
+ else if(*n<0) info = 2;
+ else if(*incx==0) info = 5;
+ if(info)
+ return xerbla_(SCALAR_SUFFIX_UP"HPR ",&info,6);
+
+ if(alpha==Scalar(0))
+ return 1;
+
+ Scalar* x_cpy = get_compact_vector(x, *n, *incx);
+
+ int code = UPLO(*uplo);
+ if(code>=2 || func[code]==0)
+ return 0;
+
+ func[code](*n, ap, x_cpy, alpha);
+
+ if(x_cpy!=x) delete[] x_cpy;
+
+ return 1;
+}
/** ZHPR2 performs the hermitian rank 2 operation
*
diff --git a/blas/sspr.f b/blas/sspr.f
deleted file mode 100644
index bae9261..0000000
--- a/blas/sspr.f
+++ /dev/null
@@ -1,202 +0,0 @@
- SUBROUTINE SSPR(UPLO,N,ALPHA,X,INCX,AP)
-* .. Scalar Arguments ..
- REAL ALPHA
- INTEGER INCX,N
- CHARACTER UPLO
-* ..
-* .. Array Arguments ..
- REAL AP(*),X(*)
-* ..
-*
-* Purpose
-* =======
-*
-* SSPR performs the symmetric rank 1 operation
-*
-* A := alpha*x*x' + A,
-*
-* where alpha is a real scalar, x is an n element vector and A is an
-* n by n symmetric matrix, supplied in packed form.
-*
-* Arguments
-* ==========
-*
-* UPLO - CHARACTER*1.
-* On entry, UPLO specifies whether the upper or lower
-* triangular part of the matrix A is supplied in the packed
-* array AP as follows:
-*
-* UPLO = 'U' or 'u' The upper triangular part of A is
-* supplied in AP.
-*
-* UPLO = 'L' or 'l' The lower triangular part of A is
-* supplied in AP.
-*
-* Unchanged on exit.
-*
-* N - INTEGER.
-* On entry, N specifies the order of the matrix A.
-* N must be at least zero.
-* Unchanged on exit.
-*
-* ALPHA - REAL .
-* On entry, ALPHA specifies the scalar alpha.
-* Unchanged on exit.
-*
-* X - REAL array of dimension at least
-* ( 1 + ( n - 1 )*abs( INCX ) ).
-* Before entry, the incremented array X must contain the n
-* element vector x.
-* Unchanged on exit.
-*
-* INCX - INTEGER.
-* On entry, INCX specifies the increment for the elements of
-* X. INCX must not be zero.
-* Unchanged on exit.
-*
-* AP - REAL array of DIMENSION at least
-* ( ( n*( n + 1 ) )/2 ).
-* Before entry with UPLO = 'U' or 'u', the array AP must
-* contain the upper triangular part of the symmetric matrix
-* packed sequentially, column by column, so that AP( 1 )
-* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
-* and a( 2, 2 ) respectively, and so on. On exit, the array
-* AP is overwritten by the upper triangular part of the
-* updated matrix.
-* Before entry with UPLO = 'L' or 'l', the array AP must
-* contain the lower triangular part of the symmetric matrix
-* packed sequentially, column by column, so that AP( 1 )
-* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
-* and a( 3, 1 ) respectively, and so on. On exit, the array
-* AP is overwritten by the lower triangular part of the
-* updated matrix.
-*
-* Further Details
-* ===============
-*
-* Level 2 Blas routine.
-*
-* -- Written on 22-October-1986.
-* Jack Dongarra, Argonne National Lab.
-* Jeremy Du Croz, Nag Central Office.
-* Sven Hammarling, Nag Central Office.
-* Richard Hanson, Sandia National Labs.
-*
-* =====================================================================
-*
-* .. Parameters ..
- REAL ZERO
- PARAMETER (ZERO=0.0E+0)
-* ..
-* .. Local Scalars ..
- REAL TEMP
- INTEGER I,INFO,IX,J,JX,K,KK,KX
-* ..
-* .. External Functions ..
- LOGICAL LSAME
- EXTERNAL LSAME
-* ..
-* .. External Subroutines ..
- EXTERNAL XERBLA
-* ..
-*
-* Test the input parameters.
-*
- INFO = 0
- IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
- INFO = 1
- ELSE IF (N.LT.0) THEN
- INFO = 2
- ELSE IF (INCX.EQ.0) THEN
- INFO = 5
- END IF
- IF (INFO.NE.0) THEN
- CALL XERBLA('SSPR ',INFO)
- RETURN
- END IF
-*
-* Quick return if possible.
-*
- IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
-*
-* Set the start point in X if the increment is not unity.
-*
- IF (INCX.LE.0) THEN
- KX = 1 - (N-1)*INCX
- ELSE IF (INCX.NE.1) THEN
- KX = 1
- END IF
-*
-* Start the operations. In this version the elements of the array AP
-* are accessed sequentially with one pass through AP.
-*
- KK = 1
- IF (LSAME(UPLO,'U')) THEN
-*
-* Form A when upper triangle is stored in AP.
-*
- IF (INCX.EQ.1) THEN
- DO 20 J = 1,N
- IF (X(J).NE.ZERO) THEN
- TEMP = ALPHA*X(J)
- K = KK
- DO 10 I = 1,J
- AP(K) = AP(K) + X(I)*TEMP
- K = K + 1
- 10 CONTINUE
- END IF
- KK = KK + J
- 20 CONTINUE
- ELSE
- JX = KX
- DO 40 J = 1,N
- IF (X(JX).NE.ZERO) THEN
- TEMP = ALPHA*X(JX)
- IX = KX
- DO 30 K = KK,KK + J - 1
- AP(K) = AP(K) + X(IX)*TEMP
- IX = IX + INCX
- 30 CONTINUE
- END IF
- JX = JX + INCX
- KK = KK + J
- 40 CONTINUE
- END IF
- ELSE
-*
-* Form A when lower triangle is stored in AP.
-*
- IF (INCX.EQ.1) THEN
- DO 60 J = 1,N
- IF (X(J).NE.ZERO) THEN
- TEMP = ALPHA*X(J)
- K = KK
- DO 50 I = J,N
- AP(K) = AP(K) + X(I)*TEMP
- K = K + 1
- 50 CONTINUE
- END IF
- KK = KK + N - J + 1
- 60 CONTINUE
- ELSE
- JX = KX
- DO 80 J = 1,N
- IF (X(JX).NE.ZERO) THEN
- TEMP = ALPHA*X(JX)
- IX = JX
- DO 70 K = KK,KK + N - J
- AP(K) = AP(K) + X(IX)*TEMP
- IX = IX + INCX
- 70 CONTINUE
- END IF
- JX = JX + INCX
- KK = KK + N - J + 1
- 80 CONTINUE
- END IF
- END IF
-*
- RETURN
-*
-* End of SSPR .
-*
- END
diff --git a/blas/zhpr.f b/blas/zhpr.f
deleted file mode 100644
index 40efbc7..0000000
--- a/blas/zhpr.f
+++ /dev/null
@@ -1,220 +0,0 @@
- SUBROUTINE ZHPR(UPLO,N,ALPHA,X,INCX,AP)
-* .. Scalar Arguments ..
- DOUBLE PRECISION ALPHA
- INTEGER INCX,N
- CHARACTER UPLO
-* ..
-* .. Array Arguments ..
- DOUBLE COMPLEX AP(*),X(*)
-* ..
-*
-* Purpose
-* =======
-*
-* ZHPR performs the hermitian rank 1 operation
-*
-* A := alpha*x*conjg( x' ) + A,
-*
-* where alpha is a real scalar, x is an n element vector and A is an
-* n by n hermitian matrix, supplied in packed form.
-*
-* Arguments
-* ==========
-*
-* UPLO - CHARACTER*1.
-* On entry, UPLO specifies whether the upper or lower
-* triangular part of the matrix A is supplied in the packed
-* array AP as follows:
-*
-* UPLO = 'U' or 'u' The upper triangular part of A is
-* supplied in AP.
-*
-* UPLO = 'L' or 'l' The lower triangular part of A is
-* supplied in AP.
-*
-* Unchanged on exit.
-*
-* N - INTEGER.
-* On entry, N specifies the order of the matrix A.
-* N must be at least zero.
-* Unchanged on exit.
-*
-* ALPHA - DOUBLE PRECISION.
-* On entry, ALPHA specifies the scalar alpha.
-* Unchanged on exit.
-*
-* X - COMPLEX*16 array of dimension at least
-* ( 1 + ( n - 1 )*abs( INCX ) ).
-* Before entry, the incremented array X must contain the n
-* element vector x.
-* Unchanged on exit.
-*
-* INCX - INTEGER.
-* On entry, INCX specifies the increment for the elements of
-* X. INCX must not be zero.
-* Unchanged on exit.
-*
-* AP - COMPLEX*16 array of DIMENSION at least
-* ( ( n*( n + 1 ) )/2 ).
-* Before entry with UPLO = 'U' or 'u', the array AP must
-* contain the upper triangular part of the hermitian matrix
-* packed sequentially, column by column, so that AP( 1 )
-* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
-* and a( 2, 2 ) respectively, and so on. On exit, the array
-* AP is overwritten by the upper triangular part of the
-* updated matrix.
-* Before entry with UPLO = 'L' or 'l', the array AP must
-* contain the lower triangular part of the hermitian matrix
-* packed sequentially, column by column, so that AP( 1 )
-* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
-* and a( 3, 1 ) respectively, and so on. On exit, the array
-* AP is overwritten by the lower triangular part of the
-* updated matrix.
-* Note that the imaginary parts of the diagonal elements need
-* not be set, they are assumed to be zero, and on exit they
-* are set to zero.
-*
-* Further Details
-* ===============
-*
-* Level 2 Blas routine.
-*
-* -- Written on 22-October-1986.
-* Jack Dongarra, Argonne National Lab.
-* Jeremy Du Croz, Nag Central Office.
-* Sven Hammarling, Nag Central Office.
-* Richard Hanson, Sandia National Labs.
-*
-* =====================================================================
-*
-* .. Parameters ..
- DOUBLE COMPLEX ZERO
- PARAMETER (ZERO= (0.0D+0,0.0D+0))
-* ..
-* .. Local Scalars ..
- DOUBLE COMPLEX TEMP
- INTEGER I,INFO,IX,J,JX,K,KK,KX
-* ..
-* .. External Functions ..
- LOGICAL LSAME
- EXTERNAL LSAME
-* ..
-* .. External Subroutines ..
- EXTERNAL XERBLA
-* ..
-* .. Intrinsic Functions ..
- INTRINSIC DBLE,DCONJG
-* ..
-*
-* Test the input parameters.
-*
- INFO = 0
- IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
- INFO = 1
- ELSE IF (N.LT.0) THEN
- INFO = 2
- ELSE IF (INCX.EQ.0) THEN
- INFO = 5
- END IF
- IF (INFO.NE.0) THEN
- CALL XERBLA('ZHPR ',INFO)
- RETURN
- END IF
-*
-* Quick return if possible.
-*
- IF ((N.EQ.0) .OR. (ALPHA.EQ.DBLE(ZERO))) RETURN
-*
-* Set the start point in X if the increment is not unity.
-*
- IF (INCX.LE.0) THEN
- KX = 1 - (N-1)*INCX
- ELSE IF (INCX.NE.1) THEN
- KX = 1
- END IF
-*
-* Start the operations. In this version the elements of the array AP
-* are accessed sequentially with one pass through AP.
-*
- KK = 1
- IF (LSAME(UPLO,'U')) THEN
-*
-* Form A when upper triangle is stored in AP.
-*
- IF (INCX.EQ.1) THEN
- DO 20 J = 1,N
- IF (X(J).NE.ZERO) THEN
- TEMP = ALPHA*DCONJG(X(J))
- K = KK
- DO 10 I = 1,J - 1
- AP(K) = AP(K) + X(I)*TEMP
- K = K + 1
- 10 CONTINUE
- AP(KK+J-1) = DBLE(AP(KK+J-1)) + DBLE(X(J)*TEMP)
- ELSE
- AP(KK+J-1) = DBLE(AP(KK+J-1))
- END IF
- KK = KK + J
- 20 CONTINUE
- ELSE
- JX = KX
- DO 40 J = 1,N
- IF (X(JX).NE.ZERO) THEN
- TEMP = ALPHA*DCONJG(X(JX))
- IX = KX
- DO 30 K = KK,KK + J - 2
- AP(K) = AP(K) + X(IX)*TEMP
- IX = IX + INCX
- 30 CONTINUE
- AP(KK+J-1) = DBLE(AP(KK+J-1)) + DBLE(X(JX)*TEMP)
- ELSE
- AP(KK+J-1) = DBLE(AP(KK+J-1))
- END IF
- JX = JX + INCX
- KK = KK + J
- 40 CONTINUE
- END IF
- ELSE
-*
-* Form A when lower triangle is stored in AP.
-*
- IF (INCX.EQ.1) THEN
- DO 60 J = 1,N
- IF (X(J).NE.ZERO) THEN
- TEMP = ALPHA*DCONJG(X(J))
- AP(KK) = DBLE(AP(KK)) + DBLE(TEMP*X(J))
- K = KK + 1
- DO 50 I = J + 1,N
- AP(K) = AP(K) + X(I)*TEMP
- K = K + 1
- 50 CONTINUE
- ELSE
- AP(KK) = DBLE(AP(KK))
- END IF
- KK = KK + N - J + 1
- 60 CONTINUE
- ELSE
- JX = KX
- DO 80 J = 1,N
- IF (X(JX).NE.ZERO) THEN
- TEMP = ALPHA*DCONJG(X(JX))
- AP(KK) = DBLE(AP(KK)) + DBLE(TEMP*X(JX))
- IX = JX
- DO 70 K = KK + 1,KK + N - J
- IX = IX + INCX
- AP(K) = AP(K) + X(IX)*TEMP
- 70 CONTINUE
- ELSE
- AP(KK) = DBLE(AP(KK))
- END IF
- JX = JX + INCX
- KK = KK + N - J + 1
- 80 CONTINUE
- END IF
- END IF
-*
- RETURN
-*
-* End of ZHPR .
-*
- END
