SUBROUTINE PDLAHRD( N, K, NB, A, IA, JA, DESCA, TAU, T, Y, IY, JY,
     $                    DESCY, WORK )
*
*  -- ScaLAPACK auxiliary routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     May 1, 1997
*
*     .. Scalar Arguments ..
      INTEGER             IA, IY, JA, JY, K, N, NB
*     ..
*     .. Array Arguments ..
      INTEGER             DESCA( * ), DESCY( * )
      DOUBLE PRECISION    A( * ), T( * ), TAU( * ), WORK( * ), Y( * )
*     ..
*
*  Purpose
*  =======
*
*  PDLAHRD reduces the first NB columns of a real general N-by-(N-K+1)
*  distributed matrix A(IA:IA+N-1,JA:JA+N-K) so that elements below the
*  k-th subdiagonal are zero. The reduction is performed by an orthogo-
*  nal similarity transformation Q' * A * Q. The routine returns the
*  matrices V and T which determine Q as a block reflector I - V*T*V',
*  and also the matrix Y = A * V * T.
*
*  This is an auxiliary routine called by PDGEHRD. In the following
*  comments sub( A ) denotes A(IA:IA+N-1,JA:JA+N-1).
*
*  Arguments
*  =========
*
*  N       (global input) INTEGER
*          The number of rows and columns to be operated on, i.e. the
*          order of the distributed submatrix sub( A ).
*          N >= 0.
*
*  K       (global input) INTEGER
*          The offset for the reduction. Elements below the k-th
*          subdiagonal in the first NB columns are reduced to zero.
*
*  NB      (global input) INTEGER
*          The number of columns to be reduced.
*
*  A       (local input/local output) DOUBLE PRECISION pointer into
*          the local memory to an array of dimension (LLD_A,
*          LOCc(JA+N-K)). On entry, this array contains the the local
*          pieces of the N-by-(N-K+1) general distributed matrix
*          A(IA:IA+N-1,JA:JA+N-K). 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 distributed
*          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(IA:IA+N-1,JA:JA+N-K) are
*          unchanged. See Further Details.
*
*  IA      (global input) INTEGER
*          The row index in the global array A indicating the first
*          row of sub( A ).
*
*  JA      (global input) INTEGER
*          The column index in the global array A indicating the
*          first column of sub( A ).
*
*  DESCA   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix A.
*
*  TAU     (local output) DOUBLE PRECISION array, dimension LOCc(JA+N-2)
*          The scalar factors of the elementary reflectors (see Further
*          Details). TAU is tied to the distributed matrix A.
*
*  T       (local output) DOUBLE PRECISION array, dimension (NB_A,NB_A)
*          The upper triangular matrix T.
*
*  Y       (local output) DOUBLE PRECISION pointer into the local memory
*          to an array of dimension (LLD_Y,NB_A). On exit, this array
*          contains the local pieces of the N-by-NB distributed
*          matrix Y. LLD_Y >= LOCr(IA+N-1).
*
*  IY      (global input) INTEGER
*          The row index in the global array Y indicating the first
*          row of sub( Y ).
*
*  JY      (global input) INTEGER
*          The column index in the global array Y indicating the
*          first column of sub( Y ).
*
*  DESCY   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix Y.
*
*  WORK    (local workspace) DOUBLE PRECISION array, dimension (NB)
*
*  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'
*
*  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(ia+i+k:ia+n-1,ja+i-1), and tau in TAU(ja+i-1).
*
*  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(ia:ia+n-1,ja:ja+n-k) := (I-V*T*V')*(A(ia:ia+n-1,ja:ja+n-k)-Y*V').
*
*  The contents of A(ia:ia+n-1,ja:ja+n-k) on exit are illustrated by the
*  following example with n = 7, k = 3 and nb = 2:
*
*     ( a   h   a   a   a )
*     ( a   h   a   a   a )
*     ( a   h   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(ia:ia+n-1,ja:ja+n-k), h denotes a modified element of the upper
*  Hessenberg matrix H, and vi denotes an element of the vector
*  defining H(i).
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
     $                   LLD_, MB_, M_, NB_, N_, RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
     $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
     $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
      DOUBLE PRECISION   ONE, ZERO
      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            IPROC
      INTEGER            I, IACOL, IAROW, ICTXT, IOFF, II, J, JJ, JL,
     $                   JT, JW, L, MYROW, MYCOL, NPCOL, NPROW, NQ
      DOUBLE PRECISION   EI
*     ..
*     .. Local Arrays ..
      INTEGER            DESCW( DLEN_ )
*     ..
*     .. External Functions ..
      INTEGER            NUMROC
      EXTERNAL           NUMROC
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, DAXPY, DESCSET, DCOPY,
     $                   DSCAL, DTRMV, INFOG2L, PDELSET,
     $                   PDGEMV, PDLARFG, PDSCAL
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MIN, MOD
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( N.LE.1 )
     $   RETURN
*
      ICTXT = DESCA( CTXT_ )
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
      IOFF = MOD( JA-1, DESCA( NB_ ) )
      CALL INFOG2L( IA+K, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, II,
     $              JJ, IAROW, IACOL )
*
      IPROC = ( MYROW.EQ.IAROW .AND. MYCOL.EQ.IACOL )
      NQ = NUMROC( N+JA-1, DESCA( NB_ ), MYCOL, IACOL, NPCOL )
      IF( MYCOL.EQ.IACOL )
     $   NQ = NQ - IOFF
*
      EI = ZERO

      JW = IOFF + 1
      CALL DESCSET( DESCW, 1, DESCA( MB_ ), 1, DESCA( MB_ ), IAROW,
     $              IACOL, ICTXT, 1 )
*
      DO 10 L = 1, NB
         I = IA + K + L - 2
         J = JA + L - 1
*
         IF( L.GT.1 ) THEN
*
*           Update A(ia:ia+n-1,j)
*
*           Compute i-th column of A - Y * V'
*
            CALL PDGEMV( 'No transpose', N, L-1, -ONE, Y, IY, JY, DESCY,
     $                   A, I, JA, DESCA, DESCA( M_ ), ONE, A, IA, J,
     $                   DESCA, 1 )
*
*           Apply I - V * T' * V' 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' * b1
*
            IF( IPROC ) THEN
               CALL DCOPY( L-1, A( (JJ+L-2)*DESCA( LLD_ )+II ), 1,
     $                     WORK( JW ), 1 )
               CALL DTRMV( 'Lower', 'Transpose', 'Unit', L-1,
     $                     A( (JJ-1)*DESCA( LLD_ )+II ), DESCA( LLD_ ),
     $                     WORK( JW ), 1 )
            END IF
*
*           w := w + V2'*b2
*
            CALL PDGEMV( 'Transpose', N-K-L+1, L-1, ONE, A, I+1, JA,
     $                   DESCA, A, I+1, J, DESCA, 1, ONE, WORK, 1, JW,
     $                   DESCW, DESCW( M_ ) )
*
*           w := T'*w
*
            IF( IPROC )
     $         CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', L-1, T,
     $                     DESCA( NB_ ), WORK( JW ), 1 )
*
*           b2 := b2 - V2*w
*
            CALL PDGEMV( 'No transpose', N-K-L+1, L-1, -ONE, A, I+1, JA,
     $                   DESCA, WORK, 1, JW, DESCW, DESCW( M_ ), ONE,
     $                   A, I+1, J, DESCA, 1 )
*
*           b1 := b1 - V1*w
*
            IF( IPROC ) THEN
               CALL DTRMV( 'Lower', 'No transpose', 'Unit', L-1,
     $                     A( (JJ-1)*DESCA( LLD_ )+II ), DESCA( LLD_ ),
     $                     WORK( JW ), 1 )
               CALL DAXPY( L-1, -ONE, WORK( JW ), 1,
     $                     A( ( JJ+L-2 )*DESCA( LLD_ )+II ), 1 )
            END IF
            CALL PDELSET( A, I, J-1, DESCA, EI )
         END IF
*
*        Generate the elementary reflector H(i) to annihilate
*        A(ia+k+i:ia+n-1,j)
*
         CALL PDLARFG( N-K-L+1, EI, I+1, J, A, MIN( I+2, N+IA-1 ), J,
     $                 DESCA, 1, TAU )
         CALL PDELSET( A, I+1, J, DESCA, ONE )
*
*        Compute  Y(iy:y+n-1,jy+l-1)
*
         CALL PDGEMV( 'No transpose', N, N-K-L+1, ONE, A, IA, J+1,
     $                DESCA, A, I+1, J, DESCA, 1, ZERO, Y, IY, JY+L-1,
     $                DESCY, 1 )
         CALL PDGEMV( 'Transpose', N-K-L+1, L-1, ONE, A, I+1, JA, DESCA,
     $                A, I+1, J, DESCA, 1, ZERO, WORK, 1, JW, DESCW,
     $                DESCW( M_ ) )
         CALL PDGEMV( 'No transpose', N, L-1, -ONE, Y, IY, JY, DESCY,
     $                WORK, 1, JW, DESCW, DESCW( M_ ), ONE, Y, IY,
     $                JY+L-1, DESCY, 1 )
         JL = MIN( JJ+L-1, JA+NQ-1 )
         CALL PDSCAL( N, TAU( JL ), Y, IY, JY+L-1, DESCY, 1 )
*
*        Compute T(1:i,i)
*
         IF( IPROC ) THEN
            JT = ( L-1 ) * DESCA( NB_ )
            CALL DSCAL( L-1, -TAU( JL ), WORK( JW ), 1 )
            CALL DCOPY( L-1, WORK( JW ), 1, T( JT+1 ), 1 )
            CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', L-1, T,
     $                  DESCA( NB_ ), T( JT+1 ), 1 )
            T( JT+L ) = TAU( JL )
         END IF
   10 CONTINUE
*
      CALL PDELSET( A, K+NB+IA-1, J, DESCA, EI )
*
      RETURN
*
*     End of PDLAHRD
*
      END