Fix rank-1 update for self-adjoint packed matrices.

7 files changed, 53 insertions, 863 deletions

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<typename Scalar, typename Index, int StorageOrder, int UpLo, bool ConjLhs, bool ConjRhs>
 struct selfadjoint_packed_rank1_update;
@@ -27,20 +21,20 @@ struct selfadjoint_packed_rank1_update;
 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 @@ struct selfadjoint_packed_rank1_update<Scalar,Index,ColMajor,UpLo,ConjLhs,ConjRh
 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 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<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 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