*  =====================================================================
*
*  Purpose
*  =======
*
*  SLAHR2 reduces the first NB columns of A real general n-BY-(n-k+1)
*  matrix A so that elements below the k-th subdiagonal are zero. The
*  reduction is performed by an orthogonal similarity transformation
*  Q**T * A * Q. The routine returns the matrices V and T which determine
*  Q as a block reflector I - V*T*V**T, and also the matrix Y = A * V * T.
*
*  This is an auxiliary routine called by SGEHRD.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the matrix A.
*
*  K       (input) INTEGER
*          The offset for the reduction. Elements below the k-th
*          subdiagonal in the first NB columns are reduced to zero.
*          K < N.
*
*  NB      (input) INTEGER
*          The number of columns to be reduced.
*
*  A       (input/output) REAL array, dimension (LDA,N-K+1)
*          On entry, the n-by-(n-k+1) general matrix A.
*          On exit, the elements on and above the k-th subdiagonal in
*          the first NB columns are overwritten with the corresponding
*          elements of the reduced matrix; the elements below the k-th
*          subdiagonal, with the array TAU, represent the matrix Q as a
*          product of elementary reflectors. The other columns of A are
*          unchanged. See Further Details.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,N).
*
*  TAU     (output) REAL array, dimension (NB)
*          The scalar factors of the elementary reflectors. See Further
*          Details.
*
*  T       (output) REAL array, dimension (LDT,NB)
*          The upper triangular matrix T.
*
*  LDT     (input) INTEGER
*          The leading dimension of the array T.  LDT >= NB.
*
*  Y       (output) REAL array, dimension (LDY,NB)
*          The n-by-nb matrix Y.
*
*  LDY     (input) INTEGER
*          The leading dimension of the array Y. LDY >= N.
*
*  Further Details
*  ===============
*
*  The matrix Q is represented as a product of nb elementary reflectors
*
*     Q = H(1) H(2) . . . H(nb).
*
*  Each H(i) has the form
*
*     H(i) = I - tau * v * v**T
*
*  where tau is a real scalar, and v is a real vector with
*  v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
*  A(i+k+1:n,i), and tau in TAU(i).
*
*  The elements of the vectors v together form the (n-k+1)-by-nb matrix
*  V which is needed, with T and Y, to apply the transformation to the
*  unreduced part of the matrix, using an update of the form:
*  A := (I - V*T*V**T) * (A - Y*V**T).
*
*  The contents of A on exit are illustrated by the following example
*  with n = 7, k = 3 and nb = 2:
*
*     ( a   a   a   a   a )
*     ( a   a   a   a   a )
*     ( a   a   a   a   a )
*     ( h   h   a   a   a )
*     ( v1  h   a   a   a )
*     ( v1  v2  a   a   a )
*     ( v1  v2  a   a   a )
*
*  where a denotes an element of the original matrix A, h denotes a
*  modified element of the upper Hessenberg matrix H, and vi denotes an
*  element of the vector defining H(i).
*
*  This subroutine is a slight modification of LAPACK-3.0's DLAHRD
*  incorporating improvements proposed by Quintana-Orti and Van de
*  Gejin. Note that the entries of A(1:K,2:NB) differ from those
*  returned by the original LAPACK-3.0's DLAHRD routine. (This
*  subroutine is not backward compatible with LAPACK-3.0's DLAHRD.)
*
*  References
*  ==========
*
*  Gregorio Quintana-Orti and Robert van de Geijn, "Improving the
*  performance of reduction to Hessenberg form," ACM Transactions on
*  Mathematical Software, 32(2):180-194, June 2006.
*
*  =====================================================================
*
      SUBROUTINE SLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
      INTEGER           I, J, K, K2, K3, LDA, LDT, LDY, N, NB
      REAL              A( LDA, * ), T( LDT, NB ), TAU( NB ),
     $                  Y( LDY, NB ), T11, TEMP, EI
      REAL              ZERO, ONE
      PARAMETER          ( ZERO = 0.0E+0, 
     $                     ONE = 1.0E+0 )
*
*     Quick return if possible
*
      IF( N.LE.1 )
     $   RETURN
*
      DO 10 I = 1, NB
         IF( I.GT.1 ) THEN
*
*           Update A(K+1:N,I)
*
*           Update I-th column of A - Y * V**T
*
*            CALL SGEMV( 'N', N-K, I-1, -ONE, Y(K+1,1), LDY,
*     $                  A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 )
            CALL SGEMV1( N-K, I-1, -ONE, Y(K+1,1), LDY,
     $                  A( K+I-1, 1 ), LDA, A( K+1, I ) )
*            DO J = 1,I-1
*               TEMP = -A(K+I-1,J)
*               DO K3 = K,N
*                  A(K3+1,I) = A(K3+1,I) + TEMP*Y(K3+1,J)
*               END DO
*            END DO
*
*           Apply I - V * T**T * V**T to this column (call it b) from the
*           left, using the last column of T as workspace
*
*           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows)
*                    ( V2 )             ( b2 )
*
*           where V1 is unit lower triangular
*
*           w := V1**T * b1
*            CALL SCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
*            CALL STRMV( 'L', 'T', 'UNIT', I-1, A( K+1, 1 ),
*     $                  LDA, T( 1, NB ), 1 )
*           w := w + V2**T * b2
*            CALL SGEMV( 'T', N-K-I+1, I-1, ONE,
*     $                  A( K+I, 1 ), LDA,
*     $                  A( K+I, I ), 1, ONE, T( 1, NB ), 1 )
*
*           w := V1**T * b1 + V2**T * b2
*
            DO 105 J = 1,I-1
               TEMP = A(K+J,I)
               DO 130 K3 = K+J+1,K+I-1
                  TEMP = TEMP + A(K3,J)*A(K3,I)
  130          CONTINUE
               DO 95 K3 = K+I,N
                  TEMP = TEMP + A(K3,J)*A(K3,I)
   95          CONTINUE
               T(J,NB) = TEMP
  105       CONTINUE
*
*           w := T**T * w
*
*            CALL STRMV( 'U', 'T', 'NON-UNIT', I-1, T, LDT,
*     $                  T( 1, NB ), 1 )
            DO 115 J = I-1,1,-1
                TEMP = T(J,NB)
                TEMP = TEMP*T(J,J)
                DO 85 K3 = J - 1,1,-1
                    TEMP = TEMP + T(K3,J)*T(K3,NB)
   85           CONTINUE
                T(J,NB) = TEMP
  115       CONTINUE
*
*           b2 := b2 - V2*w
*           b1 := b1 - V1*w
*
            DO 666 J = I-1,1,-1
              TEMP = T(J,NB)
              IF (TEMP.NE.ZERO) THEN
                 DO 56 K3 = K+I,N
                    A(K3,I) = A(K3,I) - TEMP*A(K3,J)
   56            CONTINUE
                 DO 50 K3 = I-1,J + 1,-1
                    T(K3,NB) = T(K3,NB) + TEMP*A(K3+K,J)
   50            CONTINUE
              END IF
  666       CONTINUE
            DO K3 = 1,I-1
               A(K3+K,I) = A(K3+K,I) - T(K3,NB)
            END DO
            A( K+I-1, I-1 ) = EI
         END IF
*
*        Generate the elementary reflector H(I) to annihilate
*        A(K+I+1:N,I)
*
         CALL SLARFG( N-K-I+1, A( K+I, I ), TAU( I ) )
         EI = A( K+I, I )
         A( K+I, I ) = ONE
*
*        Compute  Y(K+1:N,I)
*
         DO 155 K3 = K,N
            TEMP = ZERO
            DO 165 J = I,N-K+1
                TEMP = TEMP + A(K+J,I)*A(K3+1,J+1)
  165       CONTINUE
            Y(K3+1,I) = TEMP
  155    CONTINUE
         DO 111 J = 1,I-1
            TEMP = ZERO
            DO 191 K3 = K+I,N
               TEMP = TEMP + A(K3,J)*A(K3,I)
  191       CONTINUE
            T(J,I) = TEMP
            IF (TEMP.NE.ZERO) THEN
               DO 151 K3 = K,N
                  Y(K3+1,I) = Y(K3+1,I) - TEMP*Y(K3+1,J)
  151          CONTINUE
            END IF
  111    CONTINUE
         T11 = TAU( I )
*         CALL SSCAL( N-K, T11, Y( K+1, I ), 1 )
         DO J = K+1,N
            Y(J,I) = T11*Y(J,I)
         END DO
*
*        Compute T(1:I,I)
*
*         CALL SSCAL( I-1, -T11, T( 1, I ), 1 )
*         CALL STRMV( 'U', 'N', 'NON-UNIT', 
*     $               I-1, T, LDT, T( 1, I ), 1 )
         DO J = 1,I-1
            T(J,I) = -T11*T(J,I)
         END DO
         DO J = 1,I-1
            TEMP = T(J,I)
            DO 13 K3 = 1,J - 1
                 T(K3,I) = T(K3,I) + TEMP*T(K3,J)
   13       CONTINUE
            T(J,I) = TEMP*T(J,J)
         END DO
         T( I, I ) = TAU( I )
   10 CONTINUE
      A( K+NB, NB ) = EI
*     Compute Y(1:K,1:NB)
      DO 220 J = 1,NB
*         CALL SCOPY(K, A(1,J+1), 1, Y(1,J), 1)
         DO 201 K3 = 1,K
            Y(K3,J) = A(K3,J+1)
  201    CONTINUE
         DO 210 K2 = J + 1,NB
            TEMP = A(K2+K,J)
            IF (TEMP.NE.ZERO) THEN
               DO 200 K3 = 1,K
                  Y(K3,J) = Y(K3,J) + TEMP*A(K3,K2+1)
  200          CONTINUE
            END IF
  210    CONTINUE
  220 CONTINUE
      IF( N.GT.K+NB ) THEN
        DO 197 J = 1,NB
          DO 187 K2 = NB,N-K
            TEMP = A(K2+K+1,J)
            DO 177 K3 = 1,K
               Y(K3,J) = Y(K3,J) + TEMP*A(K3,K2+2)
  177       CONTINUE
  187     CONTINUE
  197   CONTINUE
      END IF
      DO 180 J = NB,1,-1
          TEMP = T(J,J)
          DO 150 K3 = 1,K    
              Y(K3,J) = TEMP*Y(K3,J)
  150     CONTINUE
          DO 170 K2 = 1,J - 1
              TEMP = T(K2,J)
              IF (TEMP.NE.ZERO) THEN
                  DO 160 K3 = 1,K
                      Y(K3,J) = Y(K3,J) + TEMP*Y(K3,K2)
  160             CONTINUE
              END IF
  170     CONTINUE
  180 CONTINUE
      RETURN
*     End of SLAHR2
      END

      SUBROUTINE SGEMV1(M,N,ALPHA,A,LDA,X,INCX,Y)
      INTEGER INCX,INCY,LDA,M,N
      REAL A(LDA,*),X(*),Y(*),ALPHA,TEMP
      REAL ONE,ZERO
      PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
      INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY,LENX,LENY
*
*     Start the operations. In this version the elements of A are
*     accessed sequentially with one pass through A.
*
*        Form  y := alpha*A*x + y.
      JX = 1
      DO 60 J = 1,N
         TEMP = ALPHA*X(JX)
         IF (TEMP.NE.ZERO) THEN
            DO 50 I = 1,M
               Y(I) = Y(I) + TEMP*A(I,J)
   50       CONTINUE
         END IF
         JX = JX + INCX
   60 CONTINUE
      RETURN
*     End of SGEMV .
      END
