pslaqr3.f

Go to the documentation of this file.

      RECURSIVE SUBROUTINE pslaqr3( 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( * )
      REAL               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) REAL 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) REAL 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) REAL array, dimension KBOT
*  SI      (global output) REAL 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) REAL 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) REAL 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) REAL 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) REAL 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; PSLAQR3
*          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. )
      REAL               zero, one
      PARAMETER          ( zero = 0.0, one = 1.0 )
*     ..
*     .. Local Scalars ..
      REAL               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 )
      REAL               ddum( 1 )
*     ..
*     .. External Functions ..
      REAL               slamch, pslange
      INTEGER            pilaenvx, numroc, indxg2p, iceil, blacs_pnum
      EXTERNAL           slamch, pilaenvx, numroc, indxg2p, pslange,
     $                   iceil, blacs_pnum
*     ..
*     .. External Subroutines ..
      EXTERNAL           pscopy, psgehrd, psgemm, slabad, pslacpy,
     $                   pslaqr1, slanv2, pslaqr0, pslarf, pslarfg,
     $                   pslaset, pstrord, pselget, pselset,
     $                   pslamve, blacs_gridinfo, blacs_gridmap,
     $                   blacs_gridexit, psgemr2d
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, float, 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, 'PSLAQR3', 'SV', jw, nb, -1, -1)
      par(2) = pilaenvx(ictxt, 18, 'PSLAQR3', 'SV', jw, nb, -1, -1)
      par(3) = pilaenvx(ictxt, 19, 'PSLAQR3', 'SV', jw, nb, -1, -1)
      par(4) = pilaenvx(ictxt, 20, 'PSLAQR3', 'SV', jw, nb, -1, -1)
      par(5) = pilaenvx(ictxt, 21, 'PSLAQR3', 'SV', jw, nb, -1, -1)
      par(6) = pilaenvx(ictxt, 22, 'PSLAQR3', '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 PSGEHRD and PSORMHR.
*
         taurows = numroc( 1, 1, mycol, descv(rsrc_), nprow )
         taucols = numroc( jw+iroffh, nb, mycol, descv(csrc_),
     $        npcol )
         CALL psgehrd( jw, 1, jw, t, 1, 1, desct, work, work, -1,
     $        info )
         lwk1 = int( work( 1 ) ) + taurows*taucols
*
*        Workspace query call to PSORMHR.
*
         CALL psormhr( '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 PSLAQR0.
*
         nmin = pilaenvx( ictxt, 12, 'PSLAQR3', '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 pslaqr0( .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 pslaqr1( .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 pstrord( '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 PSLARF 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.
*
         tzrows = numroc( jw+iroffh, nb, myrow, desct(rsrc_), nprow )
         tzcols = numroc( jw+iroffh, nb, mycol, desct(csrc_), npcol )
         lwk8 = 2*tzrows*tzcols
*
*        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 ) = float( lwkopt )
*
*     IWORK(1:NSEL) is used as the array SELECT for PSTRORD.
*
      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 = slamch( 'SAFE MINIMUM' )
      safmax = one / safmin
      CALL slabad( safmin, safmax )
      ulp = slamch( 'PRECISION' )
      smlnum = safmin*( float( n ) / ulp )
*
*     Setup deflation window.
*
      jw = min( nw, kbot-ktop+1 )
      kwtop = kbot - jw + 1
      IF( kwtop.EQ.ktop ) THEN
         s = zero
      ELSE
         CALL pselget( '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 pselget( 'All', '1-tree', SR( KWTOP ), H, KWTOP, KWTOP,
     $        DESCH )
         SI( KWTOP ) = ZERO
         NS = 1
         ND = 0
.LE.         IF( ABS( S )MAX( SMLNUM, ULP*ABS( SR( KWTOP ) ) ) )
     $        THEN
            NS = 0
            ND = 1
.GT.            IF( KWTOPKTOP )
     $         CALL PSELSET( H, KWTOP, KWTOP-1 , DESCH, ZERO )
         END IF
         RETURN
      END IF
*
.EQ..AND..EQ.      IF( KWTOPKTOP  KBOT-KWTOP1 ) THEN
*
*        2-by-2 deflation window: a little more to do.
*
         CALL PSELGET( 'all', '1-tree', AA, H, KWTOP, KWTOP, DESCH )
         CALL PSELGET( 'all', '1-tree', BB, H, KWTOP, KWTOP+1, DESCH )
         CALL PSELGET( 'all', '1-tree', CC, H, KWTOP+1, KWTOP, DESCH )
         CALL PSELGET( 'all', '1-tree', DD, H, KWTOP+1, KWTOP+1, DESCH )
         CALL SLANV2( AA, BB, CC, DD, SR(KWTOP), SI(KWTOP),
     $        SR(KWTOP+1), SI(KWTOP+1), CS, SN )
         NS = 0
         ND = 2
.EQ.         IF( CCZERO ) THEN
            I = KWTOP
.LE..AND.            IF( I+2N  WANTT )
     $         CALL PSROT( N-I-1, H, I, I+2, DESCH, DESCH(M_), H, I+1,
     $              I+2, DESCH, DESCH(M_), CS, SN, WORK, LWORK, INFO )
.GT.            IF( I1 )
     $         CALL PSROT( I-1, H, 1, I, DESCH, 1, H, 1, I+1, DESCH, 1,
     $              CS, SN, WORK, LWORK, INFO )
            IF( WANTZ )
     $         CALL PSROT( IHIZ-ILOZ+1, Z, ILOZ, I, DESCZ, 1, Z, ILOZ,
     $              I+1, DESCZ, 1, CS, SN, WORK, LWORK, INFO )
            CALL PSELSET( H, I, I, DESCH, AA )
            CALL PSELSET( H, I, I+1, DESCH, BB )
            CALL PSELSET( H, I+1, I, DESCH, CC )
            CALL PSELSET( H, I+1, I+1, DESCH, DD )
         END IF
         WORK( 1 ) = FLOAT( 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 PSLASET( 'all', IROFFH, JW+IROFFH, ZERO, ONE, T, 1, 1,
     $     DESCT )
      CALL PSLASET( 'all', JW, IROFFH, ZERO, ZERO, T, 1+IROFFH, 1,
     $     DESCT )
      CALL PSLACPY( 'all', 1, JW, H, KWTOP, KWTOP, DESCH, T, 1+IROFFH,
     $     1+IROFFH, DESCT )
      CALL PSLACPY( 'upper', JW-1, JW-1, H, KWTOP+1, KWTOP, DESCH, T,
     $     1+IROFFH+1, 1+IROFFH, DESCT )
.GT.      IF( JW2 )
     $   CALL PSLASET( 'lower', JW-2, JW-2, ZERO, ZERO, T, 1+IROFFH+2,
     $        1+IROFFH, DESCT )
      CALL PSLACPY( '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 PSLASET( 'all', JW+IROFFH, JW+IROFFH, ZERO, ONE, V, 1, 1,
     $     DESCV )
*
*     Compute the Schur form of T.
*
      NPMIN = PILAENVX( ICTXT, 23, 'pslaqr3', 'sv', JW, NB, NPROW,
     $     NPCOL )
      NMIN = PILAENVX( ICTXT, 12, 'pslaqr3', 'sv', JW, 1, JW, LWORK )
      NMAX = ( N-1 ) / 3
.LE..OR..GE.      IF( MIN(NPROW, NPCOL)NPMIN+1  RECLEVEL1 ) THEN
*
*        The AED window is large enough.
*        Compute the Schur decomposition with all processors.
*
.GT..AND..LE.         IF( JW+IROFFHNMIN  JW+IROFFHNMAX
.AND..LT.     $         RECLEVELRECMAX ) THEN
            CALL PSLAQR0( .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
.EQ..AND..EQ.            IF( DESCT(RSRC_)0  DESCT(CSRC_)0 ) THEN
.GT.               IF( JW+IROFFHDESCT( MB_ ) ) THEN
                  CALL PSLAQR1( .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
.EQ..AND..EQ.                  IF( MYROW0  MYCOL0 ) THEN
                     CALL SLAHQR( .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
.GT.                  IF( NPROCS1 )
     $               CALL IGAMN2D( ICTXT, 'all', '1-tree', 1, 1, INFQR,
     $                    1, -1, -1, -1, -1, -1 )
               END IF
.LE.            ELSEIF( JW+IROFFHDESCT( MB_ ) ) THEN
.EQ..AND..EQ.               IF( MYROWDESCT(RSRC_)  MYCOLDESCT(CSRC_) )
     $              THEN
                  CALL SLAHQR( .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
.GT.               IF( NPROCS1 )
     $         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 PSLAMVE( 'all', JW+IROFFH, JW+IROFFH, T, 1, 1,
     $              DESCT, WORK(IPT0), 1, 1, DESCTZ0, WORK(IPW0) )
               CALL PSLASET( 'all', JW+IROFFH, JW+IROFFH, ZERO, ONE,
     $              WORK(IPZ0), 1, 1, DESCTZ0 )
               CALL PSLAQR1( .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 PSLAMVE( 'all', JW+IROFFH, JW+IROFFH, WORK(IPT0), 1,
     $              1, DESCTZ0, T, 1, 1, DESCT, WORK(IPW0) )
               CALL PSLAMVE( '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 )
.GE..OR..GE.         IF( MYROWNPMIN  MYCOLNPMIN ) ICTXT_NEW = -1
.GE.         IF( ICTXT_NEW0 ) 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 PSGEMR2D( JW, JW, T, 1+IROFFH, 1+IROFFH, DESCT,
     $        WORK(IPT0), 1, 1, DESCTZ0, ICTXT )
.GE.         IF( ICTXT_NEW0 ) THEN
            CALL PSLASET( 'all', JW, JW, ZERO, ONE, WORK(IPZ0), 1, 1,
     $           DESCTZ0 )
            NMIN = PILAENVX( ICTXT_NEW, 12, 'pslaqr3', 'sv', JW, 1, JW,
     $           LWORK )
.GT..AND..LE..AND..LT.            IF( JWNMIN  JWNMAX  RECLEVEL1 ) THEN
               CALL PSLAQR0( .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 PSLAQR1( .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 PSGEMR2D( JW, JW, WORK(IPT0), 1, 1, DESCTZ0, T, 1+IROFFH,
     $        1+IROFFH, DESCT, ICTXT )
         CALL PSGEMR2D( JW, JW, WORK(IPZ0), 1, 1, DESCTZ0, V, 1+IROFFH,
     $        1+IROFFH, DESCV, ICTXT )
.GE.         IF( ICTXT_NEW0 )
     $      CALL BLACS_GRIDEXIT( ICTXT_NEW )
.GT.         IF( MYROW+MYCOL0 ) 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 SGSUM2D( ICTXT, 'all', ' ', JW, 1, SR(KWTOP), JW, -1, -1 )
         CALL SGSUM2D( ICTXT, 'all', ' ', JW, 1, SI(KWTOP), JW, -1, -1 )
      END IF
*
*     Adjust INFQR for offset from block border in submatrices.
*
.NE.      IF( INFQR0 )
     $   INFQR = INFQR - IROFFH
*
*     PSTRORD needs a clean margin near the diagonal.
*
      DO 50 J = 1, JW - 3
         CALL PSELSET( T, J+2, J, DESCT, ZERO )
         CALL PSELSET( T, J+3, J, DESCT, ZERO )
 50   CONTINUE
.GT.      IF( JW2 )
     $   CALL PSELSET( T, JW, JW-2, DESCT, ZERO )
*
*     Check local residual for AED Schur decomposition.
*
      RESAED = 0.0
*
*     Clean up the array SELECT for PSTRORD.
*
      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
.GT.      IF( ILST1 ) THEN
         CALL PSELGET( 'all', '1-tree', ELEM, T, ILST, ILST-1, DESCT )
.NE.         BULGE = ELEMZERO
         IF( BULGE ) ILST = ILST+1
      END IF
*
 80   CONTINUE
.LE.      IF( ILSTNS+IROFFH ) THEN
*
*        Find the top-left corner of the local window.
*
         LILST = MAX(ILST,NS+IROFFH-NB+1)
.GT.         IF( LILST1 ) THEN
            CALL PSELGET( 'all', '1-tree', ELEM, T, LILST, LILST-1,
     $           DESCT )
.NE.            BULGE = ELEMZERO
            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
.LE.         IF( LILSTNS+IROFFH ) THEN
.EQ.            IF( NS1 ) THEN
               BULGE = .FALSE.
            ELSE
               CALL PSELGET( '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 PSELGET( 'all', '1-tree', ELEM, T, NS+IROFFH,
     $              NS+IROFFH, DESCT )
               FOO = ABS( ELEM )
.EQ.               IF( FOOZERO )
     $            FOO = ABS( S )
               CALL PSELGET( '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 PSTRORD( '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 PSTRORD.
*
                  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 PSELGET( 'all', '1-tree', ELEM1, T, NS+IROFFH,
     $              NS+IROFFH, DESCT )
               CALL PSELGET( 'all', '1-tree', ELEM2, T, NS+IROFFH,
     $              NS+IROFFH-1, DESCT )
               CALL PSELGET( '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 PSELGET( 'all', '1-tree', ELEM1, V, 1+IROFFH,
     $              NS+IROFFH, DESCV )
               CALL PSELGET( '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 PSTRORD( '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 PSTRORD.
*
                  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 PSTRORD( '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 SCOPY( JW, WORK(1+IROFFH), 1, SR( KWTOP ), 1 )
      CALL SCOPY( JW, WORK(JW+2*IROFFH+1), 1, SI( KWTOP ), 1 )
*
*     Check local residual for reordered AED Schur decomposition.
*
      RESAED = 0.0
*
*     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 PSTRORD( '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 SCOPY( JW, WORK(1+IROFFH), 1, SR( KWTOP ), 1 )
               CALL SCOPY( 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 PSLASET( 'all', NS+IROFFH, 1, ZERO, ZERO, WORK(ITAU),
     $           1, 1, DESCTAU )
*
            CALL PSCOPY( NS, V, 1+IROFFH, 1+IROFFH, DESCV, DESCV(M_),
     $           WORK(IR), 1+IROFFH, 1, DESCR, 1 )
            CALL PSLARFG( NS, BETA, 1+IROFFH, 1, WORK(IR), 2+IROFFH, 1,
     $           DESCR, 1, WORK(ITAU+IROFFH) )
            CALL PSELSET( WORK(IR), 1+IROFFH, 1, DESCR, ONE )
*
            CALL PSLASET( 'lower', JW-2, JW-2, ZERO, ZERO, T, 3+IROFFH,
     $           1+IROFFH, DESCT )
*
            CALL PSLARF( 'left', NS, JW, WORK(IR), 1+IROFFH, 1, DESCR,
     $           1, WORK(ITAU+IROFFH), T, 1+IROFFH, 1+IROFFH,
     $           DESCT, WORK( IPW ) )
            CALL PSLARF( 'right', NS, NS, WORK(IR), 1+IROFFH, 1, DESCR,
     $           1, WORK(ITAU+IROFFH), T, 1+IROFFH, 1+IROFFH,
     $           DESCT, WORK( IPW ) )
            CALL PSLARF( '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 PSGEHRD( 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 PSELGET( 'all', '1-tree', ELEM, V, 1+IROFFH,
     $           1+IROFFH, DESCV )
            CALL PSELSET( H, KWTOP, KWTOP-1, DESCH, S*ELEM )
         END IF
         CALL PSLACPY( 'upper', JW-1, JW-1, T, 1+IROFFH+1, 1+IROFFH,
     $        DESCT, H, KWTOP+1, KWTOP, DESCH )
         CALL PSLACPY( 'all', 1, JW, T, 1+IROFFH, 1+IROFFH, DESCT, H,
     $        KWTOP, KWTOP, DESCH )
         CALL PSLACPY( '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 PSORMHR( '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 PSGEMM( 'no', 'No', kln, jw, jw, one, h, ltop,
     $        kwtop, desch, v, 1+iroffh, 1+iroffh, descv, zero,
     $        work, 1+iroffhh, 1+icoffhh, deschh )
         CALL pslacpy( '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 psgemm( '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 pslacpy( '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 psgemm( 'No', 'No', kln, jw, jw, one, z, iloz,
     $           kwtop, descz, v, 1+iroffh, 1+iroffh, descv,
     $           zero, work, 1+iroffzz, 1+icoffzz, desczz )
            CALL pslacpy( '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 ) = float( lwkopt )
      iwork( 1 ) = ilwkopt + nsel
*
*     End of PSLAQR3
*
      END