Comment FIXMEs on Rank2Update.h and remove unused files.

5 files changed, 19 insertions, 997 deletions

diff --git a/blas/Rank2Update.h b/blas/Rank2Update.h
index 0cf3a1961..2767916e2 100644
--- a/blas/Rank2Update.h
+++ b/blas/Rank2Update.h
@@ -23,13 +23,12 @@ struct rank2_update_selector<Scalar,Index,Upper>
 {
   static void run(Index size, Scalar* mat, Index stride, const Scalar* _u, const Scalar* _v, Scalar alpha)
   {
-    typedef Matrix<Scalar,Dynamic,1> PlainVector;
-    Map<const PlainVector> u(_u, size), v(_v, size);
-
+    Map<const Matrix<Scalar,Dynamic,1> > u(_u, size), v(_v, size);
     for (Index i=0; i<size; ++i)
     {
-      Map<PlainVector>(mat+stride*i, i+1) += conj(alpha) * conj(_u[i]) * v.head(i+1)
-					  +  alpha * conj(_v[i]) * u.head(i+1);
+      Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i, i+1) +=
+			conj(alpha) * conj(_u[i]) * v.head(i+1)
+		      + alpha * conj(_v[i]) * u.head(i+1);
     }
   }
 };
@@ -39,13 +38,12 @@ struct rank2_update_selector<Scalar,Index,Lower>
 {
   static void run(Index size, Scalar* mat, Index stride, const Scalar* _u, const Scalar* _v, Scalar alpha)
   {
-    typedef Matrix<Scalar,Dynamic,1> PlainVector;
-    Map<const PlainVector> u(_u, size), v(_v, size);
-
+    Map<const Matrix<Scalar,Dynamic,1> > u(_u, size), v(_v, size);
     for (Index i=0; i<size; ++i)
     {
-      Map<PlainVector>(mat+(stride+1)*i, size-i) += conj(alpha) * conj(_u[i]) * v.tail(size-i)
-						 +  alpha * conj(_v[i]) * u.tail(size-i);
+      Map<Matrix<Scalar,Dynamic,1> >(mat+(stride+1)*i, size-i) +=
+			conj(alpha) * conj(_u[i]) * v.tail(size-i)
+		      + alpha * conj(_v[i]) * u.tail(size-i);
     }
   }
 };
@@ -61,16 +59,16 @@ struct packed_rank2_update_selector<Scalar,Index,Upper>
 {
   static void run(Index size, Scalar* mat, const Scalar* _u, const Scalar* _v, Scalar alpha)
   {
-    typedef Matrix<Scalar,Dynamic,1> PlainVector;
-    Map<const PlainVector> u(_u, size), v(_v, size);
+    Map<const Matrix<Scalar,Dynamic,1> > u(_u, size), v(_v, size);
     Index offset = 0;
-
     for (Index i=0; i<size; ++i)
     {
-      offset += i;
-      Map<PlainVector>(mat+offset, i+1) += conj(alpha) * conj(_u[i]) * v.head(i+1)
-					+  alpha * conj(_v[i]) * u.head(i+1);
+      Map<Matrix<Scalar,Dynamic,1> >(mat+offset, i+1) +=
+			conj(alpha) * conj(_u[i]) * v.head(i+1)
+		      + alpha * conj(_v[i]) * u.head(i+1);
+      //FIXME This should be handled outside.
       mat[offset+i] = real(mat[offset+i]);
+      offset += i+1;
     }
   }
 };
@@ -80,14 +78,14 @@ struct packed_rank2_update_selector<Scalar,Index,Lower>
 {
   static void run(Index size, Scalar* mat, const Scalar* _u, const Scalar* _v, Scalar alpha)
   {
-    typedef Matrix<Scalar,Dynamic,1> PlainVector;
-    Map<const PlainVector> u(_u, size), v(_v, size);
+    Map<const Matrix<Scalar,Dynamic,1> > u(_u, size), v(_v, size);
     Index offset = 0;
-
     for (Index i=0; i<size; ++i)
     {
-      Map<PlainVector>(mat+offset, size-i) += conj(alpha) * conj(_u[i]) * v.tail(size-i)
-					   +  alpha * conj(_v[i]) * u.tail(size-i);
+      Map<Matrix<Scalar,Dynamic,1> >(mat+offset, size-i) +=
+			conj(alpha) * conj(_u[i]) * v.tail(size-i)
+		      + alpha * conj(_v[i]) * u.tail(size-i);
+      //FIXME This should be handled outside.
       mat[offset] = real(mat[offset]);
       offset += size-i;
     }
diff --git a/blas/chpr2.f b/blas/chpr2.f
deleted file mode 100644
index a0020ef3e..000000000
--- a/blas/chpr2.f
+++ /dev/null
@@ -1,255 +0,0 @@
-      SUBROUTINE CHPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
-*     .. Scalar Arguments ..
-      COMPLEX ALPHA
-      INTEGER INCX,INCY,N
-      CHARACTER UPLO
-*     ..
-*     .. Array Arguments ..
-      COMPLEX AP(*),X(*),Y(*)
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  CHPR2  performs the hermitian rank 2 operation
-*
-*     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
-*
-*  where alpha is a scalar, x and y are n element vectors 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  - COMPLEX         .
-*           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.
-*
-*  Y      - COMPLEX          array of dimension at least
-*           ( 1 + ( n - 1 )*abs( INCY ) ).
-*           Before entry, the incremented array Y must contain the n
-*           element vector y.
-*           Unchanged on exit.
-*
-*  INCY   - INTEGER.
-*           On entry, INCY specifies the increment for the elements of
-*           Y. INCY 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 TEMP1,TEMP2
-      INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
-*     ..
-*     .. 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
-      ELSE IF (INCY.EQ.0) THEN
-          INFO = 7
-      END IF
-      IF (INFO.NE.0) THEN
-          CALL XERBLA('CHPR2 ',INFO)
-          RETURN
-      END IF
-*
-*     Quick return if possible.
-*
-      IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
-*
-*     Set up the start points in X and Y if the increments are not both
-*     unity.
-*
-      IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
-          IF (INCX.GT.0) THEN
-              KX = 1
-          ELSE
-              KX = 1 - (N-1)*INCX
-          END IF
-          IF (INCY.GT.0) THEN
-              KY = 1
-          ELSE
-              KY = 1 - (N-1)*INCY
-          END IF
-          JX = KX
-          JY = KY
-      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) .AND. (INCY.EQ.1)) THEN
-              DO 20 J = 1,N
-                  IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*CONJG(Y(J))
-                      TEMP2 = CONJG(ALPHA*X(J))
-                      K = KK
-                      DO 10 I = 1,J - 1
-                          AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
-                          K = K + 1
-   10                 CONTINUE
-                      AP(KK+J-1) = REAL(AP(KK+J-1)) +
-     +                             REAL(X(J)*TEMP1+Y(J)*TEMP2)
-                  ELSE
-                      AP(KK+J-1) = REAL(AP(KK+J-1))
-                  END IF
-                  KK = KK + J
-   20         CONTINUE
-          ELSE
-              DO 40 J = 1,N
-                  IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*CONJG(Y(JY))
-                      TEMP2 = CONJG(ALPHA*X(JX))
-                      IX = KX
-                      IY = KY
-                      DO 30 K = KK,KK + J - 2
-                          AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
-                          IX = IX + INCX
-                          IY = IY + INCY
-   30                 CONTINUE
-                      AP(KK+J-1) = REAL(AP(KK+J-1)) +
-     +                             REAL(X(JX)*TEMP1+Y(JY)*TEMP2)
-                  ELSE
-                      AP(KK+J-1) = REAL(AP(KK+J-1))
-                  END IF
-                  JX = JX + INCX
-                  JY = JY + INCY
-                  KK = KK + J
-   40         CONTINUE
-          END IF
-      ELSE
-*
-*        Form  A  when lower triangle is stored in AP.
-*
-          IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
-              DO 60 J = 1,N
-                  IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*CONJG(Y(J))
-                      TEMP2 = CONJG(ALPHA*X(J))
-                      AP(KK) = REAL(AP(KK)) +
-     +                         REAL(X(J)*TEMP1+Y(J)*TEMP2)
-                      K = KK + 1
-                      DO 50 I = J + 1,N
-                          AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
-                          K = K + 1
-   50                 CONTINUE
-                  ELSE
-                      AP(KK) = REAL(AP(KK))
-                  END IF
-                  KK = KK + N - J + 1
-   60         CONTINUE
-          ELSE
-              DO 80 J = 1,N
-                  IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*CONJG(Y(JY))
-                      TEMP2 = CONJG(ALPHA*X(JX))
-                      AP(KK) = REAL(AP(KK)) +
-     +                         REAL(X(JX)*TEMP1+Y(JY)*TEMP2)
-                      IX = JX
-                      IY = JY
-                      DO 70 K = KK + 1,KK + N - J
-                          IX = IX + INCX
-                          IY = IY + INCY
-                          AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
-   70                 CONTINUE
-                  ELSE
-                      AP(KK) = REAL(AP(KK))
-                  END IF
-                  JX = JX + INCX
-                  JY = JY + INCY
-                  KK = KK + N - J + 1
-   80         CONTINUE
-          END IF
-      END IF
-*
-      RETURN
-*
-*     End of CHPR2 .
-*
-      END
diff --git a/blas/dspr2.f b/blas/dspr2.f
deleted file mode 100644
index 6f6b54a8c..000000000
--- a/blas/dspr2.f
+++ /dev/null
@@ -1,233 +0,0 @@
-      SUBROUTINE DSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
-*     .. Scalar Arguments ..
-      DOUBLE PRECISION ALPHA
-      INTEGER INCX,INCY,N
-      CHARACTER UPLO
-*     ..
-*     .. Array Arguments ..
-      DOUBLE PRECISION AP(*),X(*),Y(*)
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  DSPR2  performs the symmetric rank 2 operation
-*
-*     A := alpha*x*y' + alpha*y*x' + A,
-*
-*  where alpha is a scalar, x and y are n element vectors 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.
-*
-*  Y      - DOUBLE PRECISION array of dimension at least
-*           ( 1 + ( n - 1 )*abs( INCY ) ).
-*           Before entry, the incremented array Y must contain the n
-*           element vector y.
-*           Unchanged on exit.
-*
-*  INCY   - INTEGER.
-*           On entry, INCY specifies the increment for the elements of
-*           Y. INCY 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 TEMP1,TEMP2
-      INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
-*     ..
-*     .. 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
-      ELSE IF (INCY.EQ.0) THEN
-          INFO = 7
-      END IF
-      IF (INFO.NE.0) THEN
-          CALL XERBLA('DSPR2 ',INFO)
-          RETURN
-      END IF
-*
-*     Quick return if possible.
-*
-      IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
-*
-*     Set up the start points in X and Y if the increments are not both
-*     unity.
-*
-      IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
-          IF (INCX.GT.0) THEN
-              KX = 1
-          ELSE
-              KX = 1 - (N-1)*INCX
-          END IF
-          IF (INCY.GT.0) THEN
-              KY = 1
-          ELSE
-              KY = 1 - (N-1)*INCY
-          END IF
-          JX = KX
-          JY = KY
-      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) .AND. (INCY.EQ.1)) THEN
-              DO 20 J = 1,N
-                  IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*Y(J)
-                      TEMP2 = ALPHA*X(J)
-                      K = KK
-                      DO 10 I = 1,J
-                          AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
-                          K = K + 1
-   10                 CONTINUE
-                  END IF
-                  KK = KK + J
-   20         CONTINUE
-          ELSE
-              DO 40 J = 1,N
-                  IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*Y(JY)
-                      TEMP2 = ALPHA*X(JX)
-                      IX = KX
-                      IY = KY
-                      DO 30 K = KK,KK + J - 1
-                          AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
-                          IX = IX + INCX
-                          IY = IY + INCY
-   30                 CONTINUE
-                  END IF
-                  JX = JX + INCX
-                  JY = JY + INCY
-                  KK = KK + J
-   40         CONTINUE
-          END IF
-      ELSE
-*
-*        Form  A  when lower triangle is stored in AP.
-*
-          IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
-              DO 60 J = 1,N
-                  IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*Y(J)
-                      TEMP2 = ALPHA*X(J)
-                      K = KK
-                      DO 50 I = J,N
-                          AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
-                          K = K + 1
-   50                 CONTINUE
-                  END IF
-                  KK = KK + N - J + 1
-   60         CONTINUE
-          ELSE
-              DO 80 J = 1,N
-                  IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*Y(JY)
-                      TEMP2 = ALPHA*X(JX)
-                      IX = JX
-                      IY = JY
-                      DO 70 K = KK,KK + N - J
-                          AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
-                          IX = IX + INCX
-                          IY = IY + INCY
-   70                 CONTINUE
-                  END IF
-                  JX = JX + INCX
-                  JY = JY + INCY
-                  KK = KK + N - J + 1
-   80         CONTINUE
-          END IF
-      END IF
-*
-      RETURN
-*
-*     End of DSPR2 .
-*
-      END
diff --git a/blas/sspr2.f b/blas/sspr2.f
deleted file mode 100644
index cd27c734b..000000000
--- a/blas/sspr2.f
+++ /dev/null
@@ -1,233 +0,0 @@
-      SUBROUTINE SSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
-*     .. Scalar Arguments ..
-      REAL ALPHA
-      INTEGER INCX,INCY,N
-      CHARACTER UPLO
-*     ..
-*     .. Array Arguments ..
-      REAL AP(*),X(*),Y(*)
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  SSPR2  performs the symmetric rank 2 operation
-*
-*     A := alpha*x*y' + alpha*y*x' + A,
-*
-*  where alpha is a scalar, x and y are n element vectors 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.
-*
-*  Y      - REAL             array of dimension at least
-*           ( 1 + ( n - 1 )*abs( INCY ) ).
-*           Before entry, the incremented array Y must contain the n
-*           element vector y.
-*           Unchanged on exit.
-*
-*  INCY   - INTEGER.
-*           On entry, INCY specifies the increment for the elements of
-*           Y. INCY 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 TEMP1,TEMP2
-      INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
-*     ..
-*     .. 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
-      ELSE IF (INCY.EQ.0) THEN
-          INFO = 7
-      END IF
-      IF (INFO.NE.0) THEN
-          CALL XERBLA('SSPR2 ',INFO)
-          RETURN
-      END IF
-*
-*     Quick return if possible.
-*
-      IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
-*
-*     Set up the start points in X and Y if the increments are not both
-*     unity.
-*
-      IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
-          IF (INCX.GT.0) THEN
-              KX = 1
-          ELSE
-              KX = 1 - (N-1)*INCX
-          END IF
-          IF (INCY.GT.0) THEN
-              KY = 1
-          ELSE
-              KY = 1 - (N-1)*INCY
-          END IF
-          JX = KX
-          JY = KY
-      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) .AND. (INCY.EQ.1)) THEN
-              DO 20 J = 1,N
-                  IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*Y(J)
-                      TEMP2 = ALPHA*X(J)
-                      K = KK
-                      DO 10 I = 1,J
-                          AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
-                          K = K + 1
-   10                 CONTINUE
-                  END IF
-                  KK = KK + J
-   20         CONTINUE
-          ELSE
-              DO 40 J = 1,N
-                  IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*Y(JY)
-                      TEMP2 = ALPHA*X(JX)
-                      IX = KX
-                      IY = KY
-                      DO 30 K = KK,KK + J - 1
-                          AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
-                          IX = IX + INCX
-                          IY = IY + INCY
-   30                 CONTINUE
-                  END IF
-                  JX = JX + INCX
-                  JY = JY + INCY
-                  KK = KK + J
-   40         CONTINUE
-          END IF
-      ELSE
-*
-*        Form  A  when lower triangle is stored in AP.
-*
-          IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
-              DO 60 J = 1,N
-                  IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*Y(J)
-                      TEMP2 = ALPHA*X(J)
-                      K = KK
-                      DO 50 I = J,N
-                          AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
-                          K = K + 1
-   50                 CONTINUE
-                  END IF
-                  KK = KK + N - J + 1
-   60         CONTINUE
-          ELSE
-              DO 80 J = 1,N
-                  IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*Y(JY)
-                      TEMP2 = ALPHA*X(JX)
-                      IX = JX
-                      IY = JY
-                      DO 70 K = KK,KK + N - J
-                          AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
-                          IX = IX + INCX
-                          IY = IY + INCY
-   70                 CONTINUE
-                  END IF
-                  JX = JX + INCX
-                  JY = JY + INCY
-                  KK = KK + N - J + 1
-   80         CONTINUE
-          END IF
-      END IF
-*
-      RETURN
-*
-*     End of SSPR2 .
-*
-      END
diff --git a/blas/zhpr2.f b/blas/zhpr2.f
deleted file mode 100644
index 99977462e..000000000
--- a/blas/zhpr2.f
+++ /dev/null
@@ -1,255 +0,0 @@
-      SUBROUTINE ZHPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
-*     .. Scalar Arguments ..
-      DOUBLE COMPLEX ALPHA
-      INTEGER INCX,INCY,N
-      CHARACTER UPLO
-*     ..
-*     .. Array Arguments ..
-      DOUBLE COMPLEX AP(*),X(*),Y(*)
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  ZHPR2  performs the hermitian rank 2 operation
-*
-*     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
-*
-*  where alpha is a scalar, x and y are n element vectors 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  - COMPLEX*16      .
-*           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.
-*
-*  Y      - COMPLEX*16       array of dimension at least
-*           ( 1 + ( n - 1 )*abs( INCY ) ).
-*           Before entry, the incremented array Y must contain the n
-*           element vector y.
-*           Unchanged on exit.
-*
-*  INCY   - INTEGER.
-*           On entry, INCY specifies the increment for the elements of
-*           Y. INCY 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 TEMP1,TEMP2
-      INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
-*     ..
-*     .. 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
-      ELSE IF (INCY.EQ.0) THEN
-          INFO = 7
-      END IF
-      IF (INFO.NE.0) THEN
-          CALL XERBLA('ZHPR2 ',INFO)
-          RETURN
-      END IF
-*
-*     Quick return if possible.
-*
-      IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
-*
-*     Set up the start points in X and Y if the increments are not both
-*     unity.
-*
-      IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
-          IF (INCX.GT.0) THEN
-              KX = 1
-          ELSE
-              KX = 1 - (N-1)*INCX
-          END IF
-          IF (INCY.GT.0) THEN
-              KY = 1
-          ELSE
-              KY = 1 - (N-1)*INCY
-          END IF
-          JX = KX
-          JY = KY
-      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) .AND. (INCY.EQ.1)) THEN
-              DO 20 J = 1,N
-                  IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*DCONJG(Y(J))
-                      TEMP2 = DCONJG(ALPHA*X(J))
-                      K = KK
-                      DO 10 I = 1,J - 1
-                          AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
-                          K = K + 1
-   10                 CONTINUE
-                      AP(KK+J-1) = DBLE(AP(KK+J-1)) +
-     +                             DBLE(X(J)*TEMP1+Y(J)*TEMP2)
-                  ELSE
-                      AP(KK+J-1) = DBLE(AP(KK+J-1))
-                  END IF
-                  KK = KK + J
-   20         CONTINUE
-          ELSE
-              DO 40 J = 1,N
-                  IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*DCONJG(Y(JY))
-                      TEMP2 = DCONJG(ALPHA*X(JX))
-                      IX = KX
-                      IY = KY
-                      DO 30 K = KK,KK + J - 2
-                          AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
-                          IX = IX + INCX
-                          IY = IY + INCY
-   30                 CONTINUE
-                      AP(KK+J-1) = DBLE(AP(KK+J-1)) +
-     +                             DBLE(X(JX)*TEMP1+Y(JY)*TEMP2)
-                  ELSE
-                      AP(KK+J-1) = DBLE(AP(KK+J-1))
-                  END IF
-                  JX = JX + INCX
-                  JY = JY + INCY
-                  KK = KK + J
-   40         CONTINUE
-          END IF
-      ELSE
-*
-*        Form  A  when lower triangle is stored in AP.
-*
-          IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
-              DO 60 J = 1,N
-                  IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*DCONJG(Y(J))
-                      TEMP2 = DCONJG(ALPHA*X(J))
-                      AP(KK) = DBLE(AP(KK)) +
-     +                         DBLE(X(J)*TEMP1+Y(J)*TEMP2)
-                      K = KK + 1
-                      DO 50 I = J + 1,N
-                          AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
-                          K = K + 1
-   50                 CONTINUE
-                  ELSE
-                      AP(KK) = DBLE(AP(KK))
-                  END IF
-                  KK = KK + N - J + 1
-   60         CONTINUE
-          ELSE
-              DO 80 J = 1,N
-                  IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
-                      TEMP1 = ALPHA*DCONJG(Y(JY))
-                      TEMP2 = DCONJG(ALPHA*X(JX))
-                      AP(KK) = DBLE(AP(KK)) +
-     +                         DBLE(X(JX)*TEMP1+Y(JY)*TEMP2)
-                      IX = JX
-                      IY = JY
-                      DO 70 K = KK + 1,KK + N - J
-                          IX = IX + INCX
-                          IY = IY + INCY
-                          AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
-   70                 CONTINUE
-                  ELSE
-                      AP(KK) = DBLE(AP(KK))
-                  END IF
-                  JX = JX + INCX
-                  JY = JY + INCY
-                  KK = KK + N - J + 1
-   80         CONTINUE
-          END IF
-      END IF
-*
-      RETURN
-*
-*     End of ZHPR2 .
-*
-      END