pdlaqr3.f

Go to the documentation of this file.

      RECURSIVE SUBROUTINE pdlaqr3( WANTT, WANTZ, N, KTOP, KBOT, NW, H,
     $                              DESCH, ILOZ, IHIZ, Z, DESCZ, NS, ND,
     $                              SR, SI, V, DESCV, NH, T, DESCT, NV,
     $                              WV, DESCW, WORK, LWORK, IWORK,
     $                              LIWORK, RECLEVEL )
*
*     Contribution from the Department of Computing Science and HPC2N,
*     Umea University, Sweden
*
*  -- ScaLAPACK auxiliary routine (version 2.0.1) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     Univ. of Colorado Denver and University of California, Berkeley.
*     January, 2012
*
      IMPLICIT NONE
*
*     .. Scalar Arguments ..
      INTEGER            ihiz, ILOZ, kbot, ktop, lwork, n, nd, nh, ns,
     $                   nv, nw, liwork, reclevel
      LOGICAL            wantt, wantz
*     ..
*     .. Array Arguments ..
      INTEGER            desch( * ), descz( * ), desct( * ), descv( * ),
     $                   descw( * ), iwork( * )
      DOUBLE PRECISION   h( * ), si( kbot ), sr( kbot ), t( * ),
     $                   v( * ), work( * ), wv( * ),
     $                   z( * )
*     ..
*
*  Purpose
*  =======
*
*  Aggressive early deflation:
*
*  This subroutine accepts as input an upper Hessenberg matrix H and
*  performs an orthogonal similarity transformation designed to detect
*  and deflate fully converged eigenvalues from a trailing principal
*  submatrix.  On output H has been overwritten by a new Hessenberg
*  matrix that is a perturbation of an orthogonal similarity
*  transformation of H.  It is to be hoped that the final version of H
*  has many zero subdiagonal entries.
*
*  Notes
*  =====
*
*  Each global data object is described by an associated description
*  vector.  This vector stores the information required to establish
*  the mapping between an object element and its corresponding process
*  and memory location.
*
*  Let A be a generic term for any 2D block cyclicly distributed array.
*  Such a global array has an associated description vector DESCA.
*  In the following comments, the character _ should be read as
*  "of the global array".
*
*  NOTATION        STORED IN      EXPLANATION
*  --------------- -------------- --------------------------------------
*  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
*                                 DTYPE_A = 1.
*  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
*                                 the BLACS process grid A is distribu-
*                                 ted over. The context itself is glo-
*                                 bal, but the handle (the integer
*                                 value) may vary.
*  M_A    (global) DESCA( M_ )    The number of rows in the global
*                                 array A.
*  N_A    (global) DESCA( N_ )    The number of columns in the global
*                                 array A.
*  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
*                                 the rows of the array.
*  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
*                                 the columns of the array.
*  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
*                                 row of the array A is distributed.
*  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
*                                 first column of the array A is
*                                 distributed.
*  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
*                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
*
*  Let K be the number of rows or columns of a distributed matrix,
*  and assume that its process grid has dimension p x q.
*  LOCr( K ) denotes the number of elements of K that a process
*  would receive if K were distributed over the p processes of its
*  process column.
*  Similarly, LOCc( K ) denotes the number of elements of K that a
*  process would receive if K were distributed over the q processes of
*  its process row.
*  The values of LOCr() and LOCc() may be determined via a call to the
*  ScaLAPACK tool function, NUMROC:
*          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
*          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
*  An upper bound for these quantities may be computed by:
*          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
*          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
*
*  Arguments
*  =========
*
*  WANTT   (global input) LOGICAL
*          If .TRUE., then the Hessenberg matrix H is fully updated
*          so that the quasi-triangular Schur factor may be
*          computed (in cooperation with the calling subroutine).
*          If .FALSE., then only enough of H is updated to preserve
*          the eigenvalues.
*
*  WANTZ   (global input) LOGICAL
*          If .TRUE., then the orthogonal matrix Z is updated so
*          so that the orthogonal Schur factor may be computed
*          (in cooperation with the calling subroutine).
*          If .FALSE., then Z is not referenced.
*
*  N       (global input) INTEGER
*          The order of the matrix H and (if WANTZ is .TRUE.) the
*          order of the orthogonal matrix Z.
*
*  KTOP    (global input) INTEGER
*          It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
*          KBOT and KTOP together determine an isolated block
*          along the diagonal of the Hessenberg matrix.
*
*  KBOT    (global input) INTEGER
*          It is assumed without a check that either
*          KBOT = N or H(KBOT+1,KBOT)=0.  KBOT and KTOP together
*          determine an isolated block along the diagonal of the
*          Hessenberg matrix.
*
*  NW      (global input) INTEGER
*          Deflation window size.  1 .LE. NW .LE. (KBOT-KTOP+1).
*
*  H       (local input/output) DOUBLE PRECISION array, dimension
*             (DESCH(LLD_),*)
*          On input the initial N-by-N section of H stores the
*          Hessenberg matrix undergoing aggressive early deflation.
*          On output H has been transformed by an orthogonal
*          similarity transformation, perturbed, and the returned
*          to Hessenberg form that (it is to be hoped) has some
*          zero subdiagonal entries.
*
*  DESCH   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix H.
*
*  ILOZ    (global input) INTEGER
*  IHIZ    (global input) INTEGER
*          Specify the rows of Z to which transformations must be
*          applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
*
*  Z       (input/output) DOUBLE PRECISION array, dimension
*             (DESCH(LLD_),*)
*          IF WANTZ is .TRUE., then on output, the orthogonal
*          similarity transformation mentioned above has been
*          accumulated into Z(ILOZ:IHIZ,ILO:IHI) from the right.
*          If WANTZ is .FALSE., then Z is unreferenced.
*
*  DESCZ   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix Z.
*
*  NS      (global output) INTEGER
*          The number of unconverged (ie approximate) eigenvalues
*          returned in SR and SI that may be used as shifts by the
*          calling subroutine.
*
*  ND      (global output) INTEGER
*          The number of converged eigenvalues uncovered by this
*          subroutine.
*
*  SR      (global output) DOUBLE PRECISION array, dimension KBOT
*  SI      (global output) DOUBLE PRECISION array, dimension KBOT
*          On output, the real and imaginary parts of approximate
*          eigenvalues that may be used for shifts are stored in
*          SR(KBOT-ND-NS+1) through SR(KBOT-ND) and
*          SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively.
*          The real and imaginary parts of converged eigenvalues
*          are stored in SR(KBOT-ND+1) through SR(KBOT) and
*          SI(KBOT-ND+1) through SI(KBOT), respectively.
*
*  V       (global workspace) DOUBLE PRECISION array, dimension 
*             (DESCV(LLD_),*)
*          An NW-by-NW distributed work array.
*
*  DESCV   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix V.
*
*  NH      (input) INTEGER scalar
*          The number of columns of T.  NH.GE.NW.
*
*  T       (global workspace) DOUBLE PRECISION array, dimension 
*             (DESCV(LLD_),*)
*
*  DESCT   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix T.
*
*  NV      (global input) INTEGER
*          The number of rows of work array WV available for
*          workspace.  NV.GE.NW.
*
*  WV      (global workspace) DOUBLE PRECISION array, dimension 
*             (DESCW(LLD_),*)
*
*  DESCW   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix WV.
*
*  WORK    (local workspace) DOUBLE PRECISION array, dimension LWORK.
*          On exit, WORK(1) is set to an estimate of the optimal value
*          of LWORK for the given values of N, NW, KTOP and KBOT.
*
*  LWORK   (local input) INTEGER
*          The dimension of the work array WORK.  LWORK = 2*NW
*          suffices, but greater efficiency may result from larger
*          values of LWORK.
*
*          If LWORK = -1, then a workspace query is assumed; PDLAQR3
*          only estimates the optimal workspace size for the given
*          values of N, NW, KTOP and KBOT.  The estimate is returned
*          in WORK(1).  No error message related to LWORK is issued
*          by XERBLA.  Neither H nor Z are accessed.
*
*  IWORK   (local workspace) INTEGER array, dimension (LIWORK)
*
*  LIWORK  (local input) INTEGER
*          The length of the workspace array IWORK
*
*  ================================================================
*  Based on contributions by
*        Robert Granat and Meiyue Shao,
*        Department of Computing Science and HPC2N,
*        Umea University, Sweden
*
*  ================================================================
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, csrc_, ctxt_, dlen_, dtype_,
     $                   lld_, mb_, m_, nb_, n_, rsrc_
      INTEGER            recmax
      LOGICAL            sortgrad
      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, recmax = 3,
     $                     sortgrad = .false. )
      DOUBLE PRECISION   zero, one
      PARAMETER          ( zero = 0.0d0, one = 1.0d0 )
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION   aa, bb, beta, cc, cs, dd, evi, evk, foo, s,
     $                   safmax, safmin, smlnum, sn, tau, ulp,
     $                   elem, elem1, elem2, elem3, r1, anorm, rnorm,
     $                   resaed
      INTEGER            i, ifst, ilst, info, infqr, j, jw, k, kcol,
     $                   kend, kln, krow, kwtop, ltop, lwk1, lwk2, lwk3,
     $                   lwkopt, nmin, lldh, lldz, lldt, lldv, lldwv,
     $                   ictxt, nprow, nmax, npcol, myrow, mycol, nb,
     $                   iroffh, m, rcols, taurows, rrows, taucols,
     $                   itau, ir, ipw, nprocs, mloc, iroffhh,
     $                   icoffhh, hhrsrc, hhcsrc, hhrows, hhcols,
     $                   iroffzz, icoffzz, zzrsrc, zzcsrc, zzrows,
     $                   zzcols, ierr, tzrows0, tzcols0, ierr0, ipt0,
     $                   ipz0, ipw0, nb2, round, lilst, kk, lilst0,
     $                   iwrk1, rsrc, csrc, lwk4, lwk5, iwrk2, lwk6,
     $                   lwk7, lwk8, ilwkopt, tzrows, tzcols, nsel,
     $                   npmin, ictxt_new, myrow_new, mycol_new
      LOGICAL            bulge, sorted, lquery
*     ..
*     .. Local Arrays ..
      INTEGER            par( 6 ), descr( dlen_ ),
     $                   desctau( dlen_ ), deschh( dlen_ ),
     $                   desczz( dlen_ ), desctz0( dlen_ ),
     $                   pmap( 64*64 )
      DOUBLE PRECISION   ddum( 1 )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   dlamch, pdlange
      INTEGER            pilaenvx, numroc, indxg2p, iceil, blacs_pnum
      EXTERNAL           dlamch, pilaenvx, numroc, indxg2p, pdlange,
     $                   mpi_wtime, iceil, blacs_pnum
*     ..
*     .. External Subroutines ..
      EXTERNAL           pdcopy, pdgehrd, pdgemm, dlabad, pdlacpy,
     $                   pdlaqr1, dlanv2, pdlaqr0, pdlarf, pdlarfg,
     $                   pdlaset, pdtrord, pdelget, pdelset,
     $                   pdlamve, blacs_gridinfo, blacs_gridmap,
     $                   blacs_gridexit, pdgemr2d
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, int, max, min, sqrt
*     ..
*     .. Executable Statements ..
      ictxt = desch( ctxt_ )
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
      nprocs = nprow*npcol
*
*     Extract local leading dimensions, blockfactors, offset for
*     keeping the alignment requirements and size of deflation window.
*
      lldh  = desch( lld_ )
      lldz  = descz( lld_ )
      lldt  = desct( lld_ )
      lldv  = descv( lld_ )
      lldwv = descw( lld_ )
      nb = desch( mb_ )
      iroffh = mod( ktop - 1, nb )
      jw = min( nw, kbot-ktop+1 )
      nsel = nb+jw
*
*     Extract environment variables for parallel eigenvalue reordering.
*
      par(1) = pilaenvx(ictxt, 17, 'PDLAQR3', 'SV', jw, nb, -1, -1)
      par(2) = pilaenvx(ictxt, 18, 'PDLAQR3', 'SV', jw, nb, -1, -1)
      par(3) = pilaenvx(ictxt, 19, 'PDLAQR3', 'SV', jw, nb, -1, -1)
      par(4) = pilaenvx(ictxt, 20, 'PDLAQR3', 'SV', jw, nb, -1, -1)
      par(5) = pilaenvx(ictxt, 21, 'PDLAQR3', 'SV', jw, nb, -1, -1)
      par(6) = pilaenvx(ictxt, 22, 'PDLAQR3', 'SV', jw, nb, -1, -1)
*
*     Check if workspace query.
*
      lquery = lwork.EQ.-1 .OR. liwork.EQ.-1
*
*     Estimate optimal workspace.
*
      IF( jw.LE.2 ) THEN
         lwkopt = 1
      ELSE
*
*        Workspace query calls to PDGEHRD and PDORMHR.
*
         taurows = numroc( 1, 1, mycol, descv(rsrc_), nprow )
         taucols = numroc( jw+iroffh, nb, mycol, descv(csrc_),
     $        npcol )
         CALL pdgehrd( jw, 1, jw, t, 1, 1, desct, work, work, -1,
     $        info )
         lwk1 = int( work( 1 ) ) + taurows*taucols
*
*        Workspace query call to PDORMHR.
*
         CALL pdormhr( 'Right', 'No', jw, jw, 1, jw, t, 1, 1, desct,
     $        work, v, 1, 1, descv, work, -1, info )
         lwk2 = int( work( 1 ) )
*
*        Workspace query call to PDLAQR0.
*
         nmin = pilaenvx( ictxt, 12, 'PDLAQR3', 'SV', jw, 1, jw, lwork )
         nmax = ( n-1 ) / 3
         IF( jw+iroffh.GT.nmin .AND. jw+iroffh.LE.nmax
     $        .AND. reclevel.LT.recmax ) THEN
            CALL pdlaqr0( .true., .true., jw+iroffh, 1+iroffh,
     $           jw+iroffh, t, desct, sr, si, 1, jw, v, descv,
     $           work, -1, iwork, liwork-nsel, infqr, 
     $           reclevel+1 )
            lwk3 = int( work( 1 ) )
            iwrk1 = iwork( 1 )
         ELSE
            rsrc = desct( rsrc_ )
            csrc = desct( csrc_ )
            desct( rsrc_ ) = 0
            desct( csrc_ ) = 0
            CALL pdlaqr1( .true., .true., jw+iroffh, 1, jw+iroffh, t,
     $           desct, sr, si, 1, jw+iroffh, v, descv, work, -1,
     $           iwork, liwork-nsel, infqr )
            desct( rsrc_ ) = rsrc
            desct( csrc_ ) = csrc
            lwk3 = int( work( 1 ) )
            iwrk1 = iwork( 1 )
         END IF
*
*        Workspace in case of alignment problems.
*
         tzrows0 = numroc( jw+iroffh, nb, myrow, 0, nprow )
         tzcols0 = numroc( jw+iroffh, nb, mycol, 0, npcol )
         lwk4 = 2 * tzrows0*tzcols0
*
*        Workspace check for reordering.
*
         CALL pdtrord( 'Vectors', iwork, par, jw+iroffh, t, 1, 1,
     $        desct, v, 1, 1, descv, ddum, ddum, mloc, work, -1,
     $        iwork, liwork-nsel, info )
         lwk5 = int( work( 1 ) )
         iwrk2 = iwork( 1 )
*
*        Extra workspace for reflecting back spike
*        (workspace for PDLARF approximated for simplicity).
*
         rrows =  numroc( n+iroffh, nb, myrow, descv(rsrc_), nprow )
         rcols =  numroc( 1, 1, mycol, descv(csrc_), npcol )
         lwk6 = rrows*rcols + taurows*taucols +
     $        2*iceil(iceil(jw+iroffh,nb),nprow)*nb
     $         *iceil(iceil(jw+iroffh,nb),npcol)*nb
*
*        Extra workspace needed by PBLAS update calls
*        (also estimated for simplicity).
*
         lwk7 = max( iceil(iceil(jw,nb),nprow)*nb *
     $               iceil(iceil(n-kbot,nb),npcol)*nb,
     $               iceil(iceil(ihiz-iloz+1,nb),nprow)*nb *
     $               iceil(iceil(jw,nb),npcol)*nb,
     $               iceil(iceil(kbot-jw,nb),nprow)*nb *
     $               iceil(iceil(jw,nb),npcol)*nb )
*
*        Residual check workspace.
*
         lwk8 = 0
*
*        Optimal workspace.
*
         lwkopt = max( lwk1, lwk2, lwk3+lwk4, lwk5, lwk6, lwk7, lwk8 )
         ilwkopt = max( iwrk1, iwrk2 )
      END IF
*
*     Quick return in case of workspace query.
*
      work( 1 ) = dble( lwkopt )
*
*     IWORK(1:NSEL) is used as the array SELECT for PDTRORD.
*
      iwork( 1 ) = ilwkopt + nsel
      IF( lquery )
     $   RETURN
*
*     Nothing to do for an empty active block ...
      ns = 0
      nd = 0
      IF( ktop.GT.kbot )
     $   RETURN
*     ... nor for an empty deflation window.
*
      IF( nw.LT.1 )
     $   RETURN
*
*     Machine constants.
*
      safmin = dlamch( 'SAFE MINIMUM' )
      safmax = one / safmin
      CALL dlabad( safmin, safmax )
      ulp = dlamch( 'PRECISION' )
      smlnum = safmin*( dble( n ) / ulp )
*
*     Setup deflation window.
*
      jw = min( nw, kbot-ktop+1 )
      kwtop = kbot - jw + 1
      IF( kwtop.EQ.ktop ) THEN
         s = zero
      ELSE
         CALL pdelget( 'All', '1-Tree', s, h, kwtop, kwtop-1, desch )
      END IF
*
      IF( kbot.EQ.kwtop ) THEN
*
*        1-by-1 deflation window: not much to do.
*
         CALL pdelget( 'All', '1-Tree', sr( kwtop ), h, kwtop, kwtop,
     $        desch )
         si( kwtop ) = zero
         ns = 1
         nd = 0
         IF( abs( s ).LE.max( smlnum, ulp*abs( sr( kwtop ) ) ) )
     $        THEN
            ns = 0
            nd = 1
            IF( kwtop.GT.ktop )
     $         CALL pdelset( h, kwtop, kwtop-1 , desch, zero )
         END IF
         RETURN
      END IF
*
      IF( kwtop.EQ.ktop .AND. kbot-kwtop.EQ.1 ) THEN
*
*        2-by-2 deflation window: a little more to do.
*
         CALL pdelget( 'All', '1-Tree', aa, h, kwtop, kwtop, desch )
         CALL pdelget( 'All', '1-Tree', bb, h, kwtop, kwtop+1, desch )
         CALL pdelget( 'All', '1-Tree', cc, h, kwtop+1, kwtop, desch )
         CALL pdelget( 'All', '1-Tree', dd, h, kwtop+1, kwtop+1, desch )
         CALL dlanv2( aa, bb, cc, dd, sr(kwtop), si(kwtop),
     $        sr(kwtop+1), si(kwtop+1), cs, sn )
         ns = 0
         nd = 2
         IF( cc.EQ.zero ) THEN
            i = kwtop
            IF( i+2.LE.n .AND. wantt )
     $         CALL pdrot( n-i-1, h, i, i+2, desch, desch(m_), h, i+1,
     $              i+2, desch, desch(m_), cs, sn, work, lwork, info )
            IF( i.GT.1 )
     $         CALL pdrot( i-1, h, 1, i, desch, 1, h, 1, i+1, desch, 1,
     $              cs, sn, work, lwork, info )
            IF( wantz )
     $         CALL pdrot( ihiz-iloz+1, z, iloz, i, descz, 1, z, iloz,
     $              i+1, descz, 1, cs, sn, work, lwork, info )
            CALL pdelset( h, i, i, desch, aa )
            CALL pdelset( h, i, i+1, desch, bb )
            CALL pdelset( h, i+1, i, desch, cc )
            CALL pdelset( h, i+1, i+1, desch, dd )
         END IF
         work( 1 ) = dble( lwkopt )
         RETURN
      END IF
*
*     Calculate new value for IROFFH in case deflation window
*     was adjusted.
*
      iroffh = mod( kwtop - 1, nb )
*
*     Adjust number of rows and columns of T matrix descriptor
*     to prepare for call to PDBTRORD.
*
      desct( m_ ) = jw+iroffh
      desct( n_ ) = jw+iroffh
*
*     Convert to spike-triangular form.  (In case of a rare QR failure,
*     this routine continues to do aggressive early deflation using that
*     part of the deflation window that converged using INFQR here and
*     there to keep track.)
*
*     Copy the trailing submatrix to the working space.
*
      CALL pdlaset( 'All', iroffh, jw+iroffh, zero, one, t, 1, 1,
     $     desct )
      CALL pdlaset( 'All', jw, iroffh, zero, zero, t, 1+iroffh, 1,
     $     desct )
      CALL pdlacpy( 'All', 1, jw, h, kwtop, kwtop, desch, t, 1+iroffh,
     $     1+iroffh, desct )
      CALL pdlacpy( 'Upper', jw-1, jw-1, h, kwtop+1, kwtop, desch, t,
     $     1+iroffh+1, 1+iroffh, desct )
      IF( jw.GT.2 )
     $   CALL pdlaset( 'Lower', jw-2, jw-2, zero, zero, t, 1+iroffh+2,
     $        1+iroffh, desct )
      CALL pdlacpy( 'All', jw-1, 1, h, kwtop+1, kwtop+jw-1, desch, t,
     $     1+iroffh+1, 1+iroffh+jw-1, desct )
*
*     Initialize the working orthogonal matrix.
*
      CALL pdlaset( 'All', jw+iroffh, jw+iroffh, zero, one, v, 1, 1,
     $     descv )
*
*     Compute the Schur form of T.
*
      npmin = pilaenvx( ictxt, 23, 'PDLAQR3', 'SV', jw, nb, nprow,
     $     npcol )
      nmin = pilaenvx( ictxt, 12, 'PDLAQR3', 'SV', jw, 1, jw, lwork )
      nmax = ( n-1 ) / 3
      IF( min(nprow, npcol).LE.npmin+1 .OR. reclevel.GE.1 ) THEN
*
*        The AED window is large enough.
*        Compute the Schur decomposition with all processors.
*
         IF( jw+iroffh.GT.nmin .AND. jw+iroffh.LE.nmax
     $        .AND. reclevel.LT.recmax ) THEN
            CALL pdlaqr0( .true., .true., jw+iroffh, 1+iroffh,
     $           jw+iroffh, t, desct, sr( kwtop-iroffh ),
     $           si( kwtop-iroffh ), 1+iroffh, jw+iroffh, v, descv,
     $           work, lwork, iwork(nsel+1), liwork-nsel, infqr,
     $           reclevel+1 )
         ELSE
            IF( desct(rsrc_).EQ.0 .AND. desct(csrc_).EQ.0 ) THEN
               IF( jw+iroffh.GT.desct( mb_ ) ) THEN
                  CALL pdlaqr1( .true., .true., jw+iroffh, 1,
     $                 jw+iroffh, t, desct, sr( kwtop-iroffh ),
     $                 si( kwtop-iroffh ), 1, jw+iroffh, v,
     $                 descv, work, lwork, iwork(nsel+1), liwork-nsel,
     $                 infqr )
               ELSE
                  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
                     CALL dlahqr( .true., .true., jw+iroffh, 1+iroffh,
     $                    jw+iroffh, t, desct(lld_),
     $                    sr( kwtop-iroffh ), si( kwtop-iroffh ),
     $                    1+iroffh, jw+iroffh, v, descv(lld_), infqr )
                  ELSE
                     infqr = 0
                  END IF
                  IF( nprocs.GT.1 )
     $               CALL igamn2d( ictxt, 'All', '1-Tree', 1, 1, infqr,
     $                    1, -1, -1, -1, -1, -1 )
               END IF
            ELSEIF( jw+iroffh.LE.desct( mb_ ) ) THEN
               IF( myrow.EQ.desct(rsrc_) .AND. mycol.EQ.desct(csrc_) )
     $              THEN
                  CALL dlahqr( .true., .true., jw+iroffh, 1+iroffh,
     $                 jw+iroffh, t, desct(lld_),
     $                 sr( kwtop-iroffh ), si( kwtop-iroffh ),
     $                 1+iroffh, jw+iroffh, v, descv(lld_), infqr )
               ELSE
                  infqr = 0
               END IF
               IF( nprocs.GT.1 )
     $         CALL igamn2d( ictxt, 'All', '1-Tree', 1, 1, infqr,
     $              1, -1, -1, -1, -1, -1 )
            ELSE
               tzrows0 = numroc( jw+iroffh, nb, myrow, 0, nprow )
               tzcols0 = numroc( jw+iroffh, nb, mycol, 0, npcol )
               CALL descinit( desctz0, jw+iroffh, jw+iroffh, nb, nb, 0,
     $              0, ictxt, max(1,tzrows0), ierr0 )
               ipt0 = 1
               ipz0 = ipt0 + max(1,tzrows0)*tzcols0
               ipw0 = ipz0 + max(1,tzrows0)*tzcols0
               CALL pdlamve( 'All', jw+iroffh, jw+iroffh, t, 1, 1,
     $              desct, work(ipt0), 1, 1, desctz0, work(ipw0) )
               CALL pdlaset( 'All', jw+iroffh, jw+iroffh, zero, one,
     $              work(ipz0), 1, 1, desctz0 )
               CALL pdlaqr1( .true., .true., jw+iroffh, 1,
     $              jw+iroffh, work(ipt0), desctz0,
     $              sr( kwtop-iroffh ), si( kwtop-iroffh ),
     $              1, jw+iroffh, work(ipz0),
     $              desctz0, work(ipw0), lwork-ipw0+1, iwork(nsel+1),
     $              liwork-nsel, infqr )
               CALL pdlamve( 'All', jw+iroffh, jw+iroffh, work(ipt0), 1,
     $              1, desctz0, t, 1, 1, desct, work(ipw0) )
               CALL pdlamve( 'All', jw+iroffh, jw+iroffh, work(ipz0), 1,
     $              1, desctz0, v, 1, 1, descv, work(ipw0) )
            END IF
         END IF
      ELSE
*
*        The AED window is too small.
*        Redistribute the AED window to a subgrid
*        and do the computation on the subgrid.
*
         ictxt_new = ictxt
         DO 20 i = 0, npmin-1
            DO 10 j = 0, npmin-1
               pmap( j+1+i*npmin ) = blacs_pnum( ictxt, i, j )
 10         CONTINUE
 20      CONTINUE
         CALL blacs_gridmap( ictxt_new, pmap, npmin, npmin, npmin )
         CALL blacs_gridinfo( ictxt_new, npmin, npmin, myrow_new,
     $        mycol_new )
         IF( myrow.GE.npmin .OR. mycol.GE.npmin ) ictxt_new = -1
         IF( ictxt_new.GE.0 ) THEN
            tzrows0 = numroc( jw, nb, myrow_new, 0, npmin )
            tzcols0 = numroc( jw, nb, mycol_new, 0, npmin )
            CALL descinit( desctz0, jw, jw, nb, nb, 0,
     $           0, ictxt_new, max(1,tzrows0), ierr0 )
            ipt0 = 1
            ipz0 = ipt0 + max(1,tzrows0)*max(1,tzcols0)
            ipw0 = ipz0 + max(1,tzrows0)*max(1,tzcols0)
         ELSE
            ipt0 = 1
            ipz0 = 2
            ipw0 = 3
            desctz0( ctxt_ ) = -1
            infqr = 0
         END IF
         CALL pdgemr2d( jw, jw, t, 1+iroffh, 1+iroffh, desct,
     $        work(ipt0), 1, 1, desctz0, ictxt )
         IF( ictxt_new.GE.0 ) THEN
            CALL pdlaset( 'All', jw, jw, zero, one, work(ipz0), 1, 1,
     $           desctz0 )
            nmin = pilaenvx( ictxt_new, 12, 'PDLAQR3', 'SV', jw, 1, jw,
     $           lwork )
            IF( jw.GT.nmin .AND. jw.LE.nmax .AND. reclevel.LT.1 ) THEN
               CALL pdlaqr0( .true., .true., jw, 1, jw, work(ipt0),
     $              desctz0, sr( kwtop ), si( kwtop ), 1, jw,
     $              work(ipz0), desctz0, work(ipw0), lwork-ipw0+1,
     $              iwork(nsel+1), liwork-nsel, infqr,
     $              reclevel+1 )
            ELSE
               CALL pdlaqr1( .true., .true., jw, 1, jw, work(ipt0),
     $              desctz0, sr( kwtop ), si( kwtop ), 1, jw,
     $              work(ipz0), desctz0, work(ipw0), lwork-ipw0+1,
     $              iwork(nsel+1), liwork-nsel, infqr )
            END IF
         END IF
         CALL pdgemr2d( jw, jw, work(ipt0), 1, 1, desctz0, t, 1+iroffh,
     $        1+iroffh, desct, ictxt )
         CALL pdgemr2d( jw, jw, work(ipz0), 1, 1, desctz0, v, 1+iroffh,
     $        1+iroffh, descv, ictxt )
         IF( ictxt_new.GE.0 )
     $      CALL blacs_gridexit( ictxt_new )
         IF( myrow+mycol.GT.0 ) THEN
            DO 40 j = 0, jw-1
               sr( kwtop+j ) = zero
               si( kwtop+j ) = zero
 40         CONTINUE
         END IF
         CALL igamn2d( ictxt, 'All', '1-Tree', 1, 1, infqr, 1, -1, -1,
     $        -1, -1, -1 )
         CALL dgsum2d( ictxt, 'All', ' ', jw, 1, sr(kwtop), jw, -1, -1 )
         CALL dgsum2d( ictxt, 'All', ' ', jw, 1, si(kwtop), jw, -1, -1 )
      END IF
*
*     Adjust INFQR for offset from block border in submatrices.
*
      IF( infqr.NE.0 )
     $   infqr = infqr - iroffh
*
*     PDTRORD needs a clean margin near the diagonal.
*
      DO 50 j = 1, jw - 3
         CALL pdelset( t, j+2, j, desct, zero )
         CALL pdelset( t, j+3, j, desct, zero )
 50   CONTINUE
      IF( jw.GT.2 )
     $   CALL pdelset( t, jw, jw-2, desct, zero )
*
*     Check local residual for AED Schur decomposition.
*
      resaed = 0.0d+00
*
*     Clean up the array SELECT for PDTRORD.
*
      DO 60 j = 1, nsel
         iwork( j ) = 0
 60   CONTINUE
*
*     Set local M counter to zero.
*
      mloc = 0
*
*     Outer deflation detection loop (label 80).
*     In this loop a bunch of undeflatable eigenvalues
*     are moved simultaneously.
*
      DO 70 j = 1, iroffh + infqr
         iwork( j ) = 1
 70   CONTINUE
*
      ns = jw
      ilst = infqr + 1 + iroffh
      IF( ilst.GT.1 ) THEN
         CALL pdelget( 'All', '1-Tree', elem, t, ilst, ilst-1, desct )
         bulge = elem.NE.zero
         IF( bulge ) ilst = ilst+1
      END IF
*
 80   CONTINUE
      IF( ilst.LE.ns+iroffh ) THEN
*
*        Find the top-left corner of the local window.
*
         lilst = max(ilst,ns+iroffh-nb+1)
         IF( lilst.GT.1 ) THEN
            CALL pdelget( 'All', '1-Tree', elem, t, lilst, lilst-1,
     $           desct )
            bulge = elem.NE.zero
            IF( bulge ) lilst = lilst+1
         END IF
*
*        Lock all eigenvalues outside the local window.
*
         DO 90 j = iroffh+1, lilst-1
            iwork( j ) = 1
 90      CONTINUE
         lilst0 = lilst
*
*        Inner deflation detection loop (label 100).
*        In this loop, the undeflatable eigenvalues are moved to the
*        top-left corner of the local window.
*
 100     CONTINUE
         IF( lilst.LE.ns+iroffh ) THEN
            IF( ns.EQ.1 ) THEN
               bulge = .false.
            ELSE
               CALL pdelget( 'All', '1-tree', ELEM, T, NS+IROFFH,
     $              NS+IROFFH-1, DESCT )
.NE.               BULGE = ELEMZERO
            END IF
*
*           Small spike tip test for deflation.
*
.NOT.            IF( BULGE ) THEN
*
*              Real eigenvalue.
*
               CALL PDELGET( 'all', '1-tree', ELEM, T, NS+IROFFH,
     $              NS+IROFFH, DESCT )
               FOO = ABS( ELEM )
.EQ.               IF( FOOZERO )
     $            FOO = ABS( S )
               CALL PDELGET( 'all', '1-tree', ELEM, V, 1+IROFFH,
     $              NS+IROFFH, DESCV )
.LE.               IF( ABS( S*ELEM )MAX( SMLNUM, ULP*FOO ) ) THEN
*
*                 Deflatable.
*
                  NS = NS - 1
               ELSE
*
*                 Undeflatable: move it up out of the way.
*
                  IFST = NS
                  DO 110 J = LILST, JW+IROFFH
                     IWORK( J ) = 0
 110              CONTINUE
                  IWORK( IFST+IROFFH ) = 1
                  CALL PDTRORD( 'vectors', IWORK, PAR, JW+IROFFH, T, 1,
     $                 1, DESCT, V, 1, 1, DESCV, WORK,
     $                 WORK(JW+IROFFH+1), MLOC,
     $                 WORK(2*(JW+IROFFH)+1), LWORK-2*(JW+IROFFH),
     $                 IWORK(NSEL+1), LIWORK-NSEL, INFO )
*
*                 Adjust the array SELECT explicitly so that it does not
*                 rely on the output of PDTRORD.
*
                  IWORK( IFST+IROFFH ) = 0
                  IWORK( LILST ) = 1
                  LILST = LILST + 1
*
*                 In case of a rare exchange failure, adjust the
*                 pointers ILST and LILST to the current place to avoid
*                 unexpected behaviors.
*
.NE.                  IF( INFO0 ) THEN
                     LILST = MAX(INFO, LILST)
                     ILST = MAX(INFO, ILST)
                  END IF
               END IF
            ELSE
*
*              Complex conjugate pair.
*
               CALL PDELGET( 'all', '1-tree', ELEM1, T, NS+IROFFH,
     $              NS+IROFFH, DESCT )
               CALL PDELGET( 'all', '1-tree', ELEM2, T, NS+IROFFH,
     $              NS+IROFFH-1, DESCT )
               CALL PDELGET( 'all', '1-tree', ELEM3, T, NS+IROFFH-1,
     $              NS+IROFFH, DESCT )
               FOO = ABS( ELEM1 ) + SQRT( ABS( ELEM2 ) )*
     $              SQRT( ABS( ELEM3 ) )
.EQ.               IF( FOOZERO )
     $            FOO = ABS( S )
               CALL PDELGET( 'all', '1-tree', ELEM1, V, 1+IROFFH,
     $              NS+IROFFH, DESCV )
               CALL PDELGET( 'all', '1-tree', ELEM2, V, 1+IROFFH,
     $              NS+IROFFH-1, DESCV )
.LE.               IF( MAX( ABS( S*ELEM1 ), ABS( S*ELEM2 ) )
     $              MAX( SMLNUM, ULP*FOO ) ) THEN
*
*                 Deflatable.
*
                  NS = NS - 2
               ELSE
*
*                 Undeflatable: move them up out of the way.
*
                  IFST = NS
                  DO 120 J = LILST, JW+IROFFH
                     IWORK( J ) = 0
 120              CONTINUE
                  IWORK( IFST+IROFFH ) = 1
                  IWORK( IFST+IROFFH-1 ) = 1
                  CALL PDTRORD( 'vectors', IWORK, PAR, JW+IROFFH, T, 1,
     $                 1, DESCT, V, 1, 1, DESCV, WORK,
     $                 WORK(JW+IROFFH+1), MLOC,
     $                 WORK(2*(JW+IROFFH)+1), LWORK-2*(JW+IROFFH),
     $                 IWORK(NSEL+1), LIWORK-NSEL, INFO )
*
*                 Adjust the array SELECT explicitly so that it does not
*                 rely on the output of PDTRORD.
*
                  IWORK( IFST+IROFFH ) = 0
                  IWORK( IFST+IROFFH-1 ) = 0
                  IWORK( LILST ) = 1
                  IWORK( LILST+1 ) = 1
                  LILST = LILST + 2
*
*                 In case of a rare exchange failure, adjust the
*                 pointers ILST and LILST to the current place to avoid
*                 unexpected behaviors.
*
.NE.                  IF( INFO0 ) THEN
                     LILST = MAX(INFO, LILST)
                     ILST = MAX(INFO, ILST)
                  END IF
               END IF
            END IF
*
*           End of inner deflation detection loop.
*
            GO TO 100
         END IF
*
*        Unlock the eigenvalues outside the local window.
*        Then undeflatable eigenvalues are moved to the proper position.
*
         DO 130 J = ILST, LILST0-1
            IWORK( J ) = 0
 130     CONTINUE
         CALL PDTRORD( 'vectors', IWORK, PAR, JW+IROFFH, T, 1, 1,
     $        DESCT, V, 1, 1, DESCV, WORK, WORK(JW+IROFFH+1),
     $        M, WORK(2*(JW+IROFFH)+1), LWORK-2*(JW+IROFFH),
     $        IWORK(NSEL+1), LIWORK-NSEL, INFO )
         ILST = M + 1
*
*        In case of a rare exchange failure, adjust the pointer ILST to
*        the current place to avoid unexpected behaviors.
*
.NE.         IF( INFO0 )
     $      ILST = MAX(INFO, ILST)
*
*        End of outer deflation detection loop.
*
         GO TO 80
      END IF
 
*
*     Post-reordering step: copy output eigenvalues to output.
*
      CALL DCOPY( JW, WORK(1+IROFFH), 1, SR( KWTOP ), 1 )
      CALL DCOPY( JW, WORK(JW+2*IROFFH+1), 1, SI( KWTOP ), 1 )
*
*     Check local residual for reordered AED Schur decomposition.
*
      RESAED = 0.0D+00
*
*     Return to Hessenberg form.
*
.EQ.      IF( NS0 )
     $   S = ZERO
*
.LT..AND.      IF( NSJW  SORTGRAD ) THEN
*
*        Sorting diagonal blocks of T improves accuracy for
*        graded matrices.  Bubble sort deals well with exchange
*        failures. Eigenvalues/shifts from T are also restored.
*
         ROUND = 0
         SORTED = .FALSE.
         I = NS + 1 + IROFFH
 140     CONTINUE
         IF( SORTED )
     $      GO TO 180
         SORTED = .TRUE.
         ROUND = ROUND + 1
*
         KEND = I - 1
         I = INFQR + 1 + IROFFH
.EQ.         IF( INS+IROFFH ) THEN
            K = I + 1
.EQ.         ELSE IF( SI( KWTOP-IROFFH + I-1 )ZERO ) THEN
            K = I + 1
         ELSE
            K = I + 2
         END IF
 150     CONTINUE
.LE.         IF( KKEND ) THEN
.EQ.            IF( KI+1 ) THEN
               EVI = ABS( SR( KWTOP-IROFFH+I-1 ) )
            ELSE
               EVI = ABS( SR( KWTOP-IROFFH+I-1 ) ) +
     $              ABS( SI( KWTOP-IROFFH+I-1 ) )
            END IF
*
.EQ.            IF( KKEND ) THEN
               EVK = ABS( SR( KWTOP-IROFFH+K-1 ) )
.EQ.            ELSEIF( SI( KWTOP-IROFFH+K-1 )ZERO ) THEN
               EVK = ABS( SR( KWTOP-IROFFH+K-1 ) )
            ELSE
               EVK = ABS( SR( KWTOP-IROFFH+K-1 ) ) +
     $              ABS( SI( KWTOP-IROFFH+K-1 ) )
            END IF
*
.GE.            IF( EVIEVK ) THEN
               I = K
            ELSE
               MLOC = 0
               SORTED = .FALSE.
               IFST = I
               ILST = K
               DO 160 J = 1, I-1
                  IWORK( J ) = 1
                  MLOC = MLOC + 1
 160           CONTINUE
.EQ.               IF( KI+2 ) THEN
                  IWORK( I ) = 0
                  IWORK(I+1) = 0
               ELSE
                  IWORK( I ) = 0
               END IF
.NE..AND..NE.               IF( KKEND  SI( KWTOP-IROFFH+K-1 )ZERO ) THEN
                  IWORK( K ) = 1
                  IWORK(K+1) = 1
                  MLOC = MLOC + 2
               ELSE
                  IWORK( K ) = 1
.LT.                  IF( KKEND ) IWORK(K+1) = 0
                  MLOC = MLOC + 1
               END IF
               DO 170 J = K+2, JW+IROFFH
                  IWORK( J ) = 0
 170           CONTINUE
               CALL PDTRORD( 'vectors', IWORK, PAR, JW+IROFFH, T, 1, 1,
     $              DESCT, V, 1, 1, DESCV, WORK, WORK(JW+IROFFH+1), M,
     $              WORK(2*(JW+IROFFH)+1), LWORK-2*(JW+IROFFH),
     $              IWORK(NSEL+1), LIWORK-NSEL, IERR )
               CALL DCOPY( JW, WORK(1+IROFFH), 1, SR( KWTOP ), 1 )
               CALL DCOPY( JW, WORK(JW+2*IROFFH+1), 1, SI( KWTOP ), 1 )
.EQ.               IF( IERR0 ) THEN
                  I = ILST
               ELSE
                  I = K
               END IF
            END IF
.EQ.            IF( IKEND ) THEN
               K = I + 1
.EQ.            ELSE IF( SI( KWTOP-IROFFH+I-1 )ZERO ) THEN
               K = I + 1
            ELSE
               K = I + 2
            END IF
            GO TO 150
         END IF
         GO TO 140
 180     CONTINUE
      END IF
*
*     Restore number of rows and columns of T matrix descriptor.
*
      DESCT( M_ ) = NW+IROFFH
      DESCT( N_ ) = NH+IROFFH
*
.LT..OR..EQ.      IF( NSJW  SZERO ) THEN
.GT..AND..NE.         IF( NS1  SZERO ) THEN
*
*           Reflect spike back into lower triangle.
*
            RROWS = NUMROC( NS+IROFFH, NB, MYROW, DESCV(RSRC_), NPROW )
            RCOLS = NUMROC( 1, 1, MYCOL, DESCV(CSRC_), NPCOL )
            CALL DESCINIT( DESCR, NS+IROFFH, 1, NB, 1, DESCV(RSRC_),
     $           DESCV(CSRC_), ICTXT, MAX(1, RROWS), INFO )
            TAUROWS = NUMROC( 1, 1, MYCOL, DESCV(RSRC_), NPROW )
            TAUCOLS = NUMROC( JW+IROFFH, NB, MYCOL, DESCV(CSRC_),
     $           NPCOL )
            CALL DESCINIT( DESCTAU, 1, JW+IROFFH, 1, NB, DESCV(RSRC_),
     $           DESCV(CSRC_), ICTXT, MAX(1, TAUROWS), INFO )
*
            IR = 1
            ITAU = IR + DESCR( LLD_ ) * RCOLS
            IPW  = ITAU + DESCTAU( LLD_ ) * TAUCOLS
*
            CALL PDLASET( 'all', NS+IROFFH, 1, ZERO, ZERO, WORK(ITAU),
     $           1, 1, DESCTAU )
*
            CALL PDCOPY( NS, V, 1+IROFFH, 1+IROFFH, DESCV, DESCV(M_),
     $           WORK(IR), 1+IROFFH, 1, DESCR, 1 )
            CALL PDLARFG( NS, BETA, 1+IROFFH, 1, WORK(IR), 2+IROFFH, 1,
     $           DESCR, 1, WORK(ITAU+IROFFH) )
            CALL PDELSET( WORK(IR), 1+IROFFH, 1, DESCR, ONE )
*
            CALL PDLASET( 'lower', JW-2, JW-2, ZERO, ZERO, T, 3+IROFFH,
     $           1+IROFFH, DESCT )
*
            CALL PDLARF( 'left', NS, JW, WORK(IR), 1+IROFFH, 1, DESCR,
     $           1, WORK(ITAU+IROFFH), T, 1+IROFFH, 1+IROFFH,
     $           DESCT, WORK( IPW ) )
            CALL PDLARF( 'right', NS, NS, WORK(IR), 1+IROFFH, 1, DESCR,
     $           1, WORK(ITAU+IROFFH), T, 1+IROFFH, 1+IROFFH,
     $           DESCT, WORK( IPW ) )
            CALL PDLARF( 'right', JW, NS, WORK(IR), 1+IROFFH, 1, DESCR,
     $           1, WORK(ITAU+IROFFH), V, 1+IROFFH, 1+IROFFH,
     $           DESCV, WORK( IPW ) )
*
            ITAU = 1
            IPW = ITAU + DESCTAU( LLD_ ) * TAUCOLS
            CALL PDGEHRD( JW+IROFFH, 1+IROFFH, NS+IROFFH, T, 1, 1,
     $           DESCT, WORK(ITAU), WORK( IPW ), LWORK-IPW+1, INFO )
         END IF
*
*        Copy updated reduced window into place.
*
.GT.         IF( KWTOP1 ) THEN
            CALL PDELGET( 'all', '1-tree', ELEM, V, 1+IROFFH,
     $           1+IROFFH, DESCV )
            CALL PDELSET( H, KWTOP, KWTOP-1, DESCH, S*ELEM )
         END IF
         CALL PDLACPY( 'upper', JW-1, JW-1, T, 1+IROFFH+1, 1+IROFFH,
     $        DESCT, H, KWTOP+1, KWTOP, DESCH )
         CALL PDLACPY( 'all', 1, JW, T, 1+IROFFH, 1+IROFFH, DESCT, H,
     $        KWTOP, KWTOP, DESCH )
         CALL PDLACPY( 'all', JW-1, 1, T, 1+IROFFH+1, 1+IROFFH+JW-1,
     $        DESCT, H, KWTOP+1, KWTOP+JW-1, DESCH )
*
*        Accumulate orthogonal matrix in order to update
*        H and Z, if requested.
*
.GT..AND..NE.         IF( NS1  SZERO ) THEN
            CALL PDORMHR( 'right', 'no', JW+IROFFH, NS+IROFFH, 1+IROFFH,
     $           NS+IROFFH, T, 1, 1, DESCT, WORK(ITAU), V, 1,
     $           1, DESCV, WORK( IPW ), LWORK-IPW+1, INFO )
         END IF
*
*        Update vertical slab in H.
*
         IF( WANTT ) THEN
            LTOP = 1
         ELSE
            LTOP = KTOP
         END IF
         KLN = MAX( 0, KWTOP-LTOP )
         IROFFHH = MOD( LTOP-1, NB )
         ICOFFHH = MOD( KWTOP-1, NB )
         HHRSRC = INDXG2P( LTOP, NB, MYROW, DESCH(RSRC_), NPROW )
         HHCSRC = INDXG2P( KWTOP, NB, MYCOL, DESCH(CSRC_), NPCOL )
         HHROWS = NUMROC( KLN+IROFFHH, NB, MYROW, HHRSRC, NPROW )
         HHCOLS = NUMROC( JW+ICOFFHH, NB, MYCOL, HHCSRC, NPCOL )
         CALL DESCINIT( DESCHH, KLN+IROFFHH, JW+ICOFFHH, NB, NB,
     $        HHRSRC, HHCSRC, ICTXT, MAX(1, HHROWS), IERR )
         CALL PDGEMM( 'no', 'no', KLN, JW, JW, ONE, H, LTOP,
     $        KWTOP, DESCH, V, 1+IROFFH, 1+IROFFH, DESCV, ZERO,
     $        WORK, 1+IROFFHH, 1+ICOFFHH, DESCHH )
         CALL PDLACPY( 'all', KLN, JW, WORK, 1+IROFFHH, 1+ICOFFHH,
     $        DESCHH, H, LTOP, KWTOP, DESCH )
*
*        Update horizontal slab in H.
*
         IF( WANTT ) THEN
            KLN = N-KBOT
            IROFFHH = MOD( KWTOP-1, NB )
            ICOFFHH = MOD( KBOT, NB )
            HHRSRC = INDXG2P( KWTOP, NB, MYROW, DESCH(RSRC_), NPROW )
            HHCSRC = INDXG2P( KBOT+1, NB, MYCOL, DESCH(CSRC_), NPCOL )
            HHROWS = NUMROC( JW+IROFFHH, NB, MYROW, HHRSRC, NPROW )
            HHCOLS = NUMROC( KLN+ICOFFHH, NB, MYCOL, HHCSRC, NPCOL )
            CALL DESCINIT( DESCHH, JW+IROFFHH, KLN+ICOFFHH, NB, NB,
     $           HHRSRC, HHCSRC, ICTXT, MAX(1, HHROWS), IERR )
            CALL PDGEMM( 'tr', 'no', JW, KLN, JW, ONE, V,
     $           1+IROFFH, 1+IROFFH, DESCV, H, KWTOP, KBOT+1,
     $           DESCH, ZERO, WORK, 1+IROFFHH, 1+ICOFFHH, DESCHH )
            CALL PDLACPY( 'all', JW, KLN, WORK, 1+IROFFHH, 1+ICOFFHH,
     $           DESCHH, H, KWTOP, KBOT+1, DESCH )
         END IF
*
*        Update vertical slab in Z.
*
         IF( WANTZ ) THEN
            KLN = IHIZ-ILOZ+1
            IROFFZZ = MOD( ILOZ-1, NB )
            ICOFFZZ = MOD( KWTOP-1, NB )
            ZZRSRC = INDXG2P( ILOZ, NB, MYROW, DESCZ(RSRC_), NPROW )
            ZZCSRC = INDXG2P( KWTOP, NB, MYCOL, DESCZ(CSRC_), NPCOL )
            ZZROWS = NUMROC( KLN+IROFFZZ, NB, MYROW, ZZRSRC, NPROW )
            ZZCOLS = NUMROC( JW+ICOFFZZ, NB, MYCOL, ZZCSRC, NPCOL )
            CALL DESCINIT( DESCZZ, KLN+IROFFZZ, JW+ICOFFZZ, NB, NB,
     $           ZZRSRC, ZZCSRC, ICTXT, MAX(1, ZZROWS), IERR )
            CALL PDGEMM( 'no', 'no', KLN, JW, JW, ONE, Z, ILOZ,
     $           KWTOP, DESCZ, V, 1+IROFFH, 1+IROFFH, DESCV,
     $           ZERO, WORK, 1+IROFFZZ, 1+ICOFFZZ, DESCZZ )
            CALL PDLACPY( 'all', KLN, JW, WORK, 1+IROFFZZ, 1+ICOFFZZ,
     $           DESCZZ, Z, ILOZ, KWTOP, DESCZ )
         END IF
      END IF
*
*     Return the number of deflations (ND) and the number of shifts (NS).
*     (Subtracting INFQR from the spike length takes care of the case of
*     a rare QR failure while calculating eigenvalues of the deflation
*     window.)
*
      ND = JW - NS
      NS = NS - INFQR
*
*     Return optimal workspace.
*
      WORK( 1 ) = DBLE( LWKOPT )
      IWORK( 1 ) = ILWKOPT + NSEL
*
*     End of PDLAQR3
*
      END