From 04f315d692d4b2b01635981d34b44743ba63571c Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Mon, 10 Sep 2012 18:25:30 +0800 Subject: Fix rank-1 update for self-adjoint packed matrices. --- blas/CMakeLists.txt | 6 +- blas/PackedSelfadjointProduct.h | 19 ++-- blas/chpr.f | 220 ---------------------------------------- blas/dspr.f | 202 ------------------------------------ blas/level2_cplx_impl.h | 47 ++++++++- blas/sspr.f | 202 ------------------------------------ blas/zhpr.f | 220 ---------------------------------------- 7 files changed, 53 insertions(+), 863 deletions(-) delete mode 100644 blas/chpr.f delete mode 100644 blas/dspr.f delete mode 100644 blas/sspr.f delete mode 100644 blas/zhpr.f (limited to 'blas') diff --git a/blas/CMakeLists.txt b/blas/CMakeLists.txt index 3877e1285..c35a2fdbe 100644 --- a/blas/CMakeLists.txt +++ b/blas/CMakeLists.txt @@ -18,9 +18,9 @@ if(EIGEN_Fortran_COMPILER_WORKS) 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 1ba67b9c1..f7c9b9341 100644 --- a/blas/PackedSelfadjointProduct.h +++ b/blas/PackedSelfadjointProduct.h @@ -14,12 +14,6 @@ namespace internal { /* 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 struct selfadjoint_packed_rank1_update; @@ -27,20 +21,20 @@ struct selfadjoint_packed_rank1_update; template struct selfadjoint_packed_rank1_update { - static void run(Index size, Scalar* mat, const Scalar* vec, Scalar alpha) + typedef typename NumTraits::Real RealScalar; + static void run(Index size, Scalar* mat, const Scalar* vec, RealScalar alpha) { typedef Map > OtherMap; typedef typename conj_expr_if::type ConjRhsType; conj_if cj; - Index offset = 0; for (Index i=0; i >(mat+offset, UpLo==Lower ? size-i : (i+1)) + Map >(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 @@ struct selfadjoint_packed_rank1_update struct selfadjoint_packed_rank1_update { - static void run(Index size, Scalar* mat, const Scalar* vec, Scalar alpha) + typedef typename NumTraits::Real RealScalar; + static void run(Index size, Scalar* mat, const Scalar* vec, RealScalar alpha) { selfadjoint_packed_rank1_update::run(size,mat,vec,alpha); } diff --git a/blas/chpr.f b/blas/chpr.f deleted file mode 100644 index 11bd5c6ee..000000000 --- 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 538e4f76b..000000000 --- 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 11ee13b4c..f52d384a9 100644 --- a/blas/level2_cplx_impl.h +++ b/blas/level2_cplx_impl.h @@ -108,10 +108,49 @@ int EIGEN_BLAS_FUNC(hemv)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa * 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::run); + func[LO] = (internal::selfadjoint_packed_rank1_update::run); + + init = true; + } + + Scalar* x = reinterpret_cast(px); + Scalar* ap = reinterpret_cast(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 bae92612e..000000000 --- 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 40efbc7d5..000000000 --- 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 -- cgit v1.2.3