pdhseqr.f

Go to the documentation of this file.

      SUBROUTINE pdhseqr( JOB, COMPZ, N, ILO, IHI, H, DESCH, WR, WI, Z,
     $                    DESCZ, WORK, LWORK, IWORK, LIWORK, INFO )
*
*     Contribution from the Department of Computing Science and HPC2N,
*     Umea University, Sweden
*
*  -- ScaLAPACK driver 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            IHI, ILO, INFO, LWORK, LIWORK, N
      CHARACTER          COMPZ, JOB
*     ..
*     .. Array Arguments ..
      INTEGER            DESCH( * ) , DESCZ( * ), IWORK( * )
      DOUBLE PRECISION   H( * ), WI( N ), WORK( * ), WR( N ), Z( * )
*     ..
*  Purpose
*  =======
*
*  PDHSEQR computes the eigenvalues of an upper Hessenberg matrix H
*  and, optionally, the matrices T and Z from the Schur decomposition
*  H = Z*T*Z**T, where T is an upper quasi-triangular matrix (the
*  Schur form), and Z is the orthogonal matrix of Schur vectors.
*
*  Optionally Z may be postmultiplied into an input orthogonal
*  matrix Q so that this routine can give the Schur factorization
*  of a matrix A which has been reduced to the Hessenberg form H
*  by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
*
*  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
*  =========
*
*  JOB     (global input) CHARACTER*1
*          = 'E':  compute eigenvalues only;
*          = 'S':  compute eigenvalues and the Schur form T.
*
*  COMPZ   (global input) CHARACTER*1
*          = 'N':  no Schur vectors are computed;
*          = 'I':  Z is initialized to the unit matrix and the matrix Z
*                  of Schur vectors of H is returned;
*          = 'V':  Z must contain an orthogonal matrix Q on entry, and
*                  the product Q*Z is returned.
*
*  N       (global input) INTEGER
*          The order of the Hessenberg matrix H (and Z if WANTZ).
*          N >= 0.
*
*  ILO     (global input) INTEGER
*  IHI     (global input) INTEGER
*          It is assumed that H is already upper triangular in rows
*          and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
*          set by a previous call to PDGEBAL, and then passed to PDGEHRD
*          when the matrix output by PDGEBAL is reduced to Hessenberg
*          form. Otherwise ILO and IHI should be set to 1 and N
*          respectively.  If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
*          If N = 0, then ILO = 1 and IHI = 0.
*
*  H       (global input/output) DOUBLE PRECISION array, dimension
*          (DESCH(LLD_),*)
*          On entry, the upper Hessenberg matrix H.
*          On exit, if JOB = 'S', H is upper quasi-triangular in
*          rows and columns ILO:IHI, with 1-by-1 and 2-by-2 blocks on
*          the main diagonal.  The 2-by-2 diagonal blocks (corresponding
*          to complex conjugate pairs of eigenvalues) are returned in
*          standard form, with H(i,i) = H(i+1,i+1) and
*          H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the
*          contents of H are unspecified on exit.
*
*  DESCH   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix H.
*
*  WR      (global output) DOUBLE PRECISION array, dimension (N)
*  WI      (global output) DOUBLE PRECISION array, dimension (N)
*          The real and imaginary parts, respectively, of the computed
*          eigenvalues ILO to IHI are stored in the corresponding
*          elements of WR and WI. If two eigenvalues are computed as a
*          complex conjugate pair, they are stored in consecutive
*          elements of WR and WI, say the i-th and (i+1)th, with
*          WI(i) > 0 and WI(i+1) < 0. If JOB = 'S', the
*          eigenvalues are stored in the same order as on the diagonal
*          of the Schur form returned in H.
*
*  Z       (global input/output) DOUBLE PRECISION array.
*          If COMPZ = 'V', on entry Z must contain the current
*          matrix Z of accumulated transformations from, e.g., PDGEHRD,
*          and on exit Z has been updated; transformations are applied
*          only to the submatrix Z(ILO:IHI,ILO:IHI).
*          If COMPZ = 'N', Z is not referenced.
*          If COMPZ = 'I', on entry Z need not be set and on exit,
*          if INFO = 0, Z contains the orthogonal matrix Z of the Schur
*          vectors of H.
*
*  DESCZ   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix Z.
*
*  WORK    (local workspace) DOUBLE PRECISION array, dimension(LWORK)
*
*  LWORK   (local input) INTEGER
*          The length of the workspace array WORK.
*
*  IWORK   (local workspace) INTEGER array, dimension (LIWORK)
*
*  LIWORK  (local input) INTEGER
*          The length of the workspace array IWORK.
*
*  INFO    (output) INTEGER
*          =    0:  successful exit
*          .LT. 0:  if INFO = -i, the i-th argument had an illegal
*                   value (see also below for -7777 and -8888).
*          .GT. 0:  if INFO = i, PDHSEQR failed to compute all of
*                   the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
*                   and WI contain those eigenvalues which have been
*                   successfully computed.  (Failures are rare.)
*
*                If INFO .GT. 0 and JOB = 'E', then on exit, the
*                remaining unconverged eigenvalues are the eigen-
*                values of the upper Hessenberg matrix rows and
*                columns ILO through INFO of the final, output
*                value of H.
*
*                If INFO .GT. 0 and JOB   = 'S', then on exit
*
*           (*)  (initial value of H)*U  = U*(final value of H)
*
*                where U is an orthogonal matrix.  The final
*                value of H is upper Hessenberg and quasi-triangular
*                in rows and columns INFO+1 through IHI.
*
*                If INFO .GT. 0 and COMPZ = 'V', then on exit
*
*                  (final value of Z)  =  (initial value of Z)*U
*
*                where U is the orthogonal matrix in (*) (regard-
*                less of the value of JOB.)
*
*                If INFO .GT. 0 and COMPZ = 'I', then on exit
*                      (final value of Z)  = U
*                where U is the orthogonal matrix in (*) (regard-
*                less of the value of JOB.)
*
*                If INFO .GT. 0 and COMPZ = 'N', then Z is not
*                accessed.
*
*          = -7777: PDLAQR0 failed to converge and PDLAQR1 was called
*                   instead. This could happen. Mostly due to a bug.
*                   Please, send a bug report to the authors.
*          = -8888: PDLAQR1 failed to converge and PDLAQR0 was called
*                   instead. This should not happen.
*
*     ================================================================
*     Based on contributions by
*        Robert Granat, Department of Computing Science and HPC2N,
*        Umea University, Sweden.
*     ================================================================
*
*     Restrictions: The block size in H and Z must be square and larger
*     than or equal to six (6) due to restrictions in PDLAQR1, PDLAQR5
*     and DLAQR6. Moreover, H and Z need to be distributed identically
*     with the same context.
*
*     ================================================================
*     References:
*       K. Braman, R. Byers, and R. Mathias,
*       The Multi-Shift QR Algorithm Part I: Maintaining Well Focused
*       Shifts, and Level 3 Performance.
*       SIAM J. Matrix Anal. Appl., 23(4):929--947, 2002.
*
*       K. Braman, R. Byers, and R. Mathias,
*       The Multi-Shift QR Algorithm Part II: Aggressive Early
*       Deflation.
*       SIAM J. Matrix Anal. Appl., 23(4):948--973, 2002.
*
*       R. Granat, B. Kagstrom, and D. Kressner,
*       A Novel Parallel QR Algorithm for Hybrid Distributed Momory HPC
*       Systems.
*       SIAM J. Sci. Comput., 32(4):2345--2378, 2010.
*
*     ================================================================
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
     $                   lld_, mb_, m_, nb_, n_, rsrc_
      LOGICAL            CRSOVER
      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,
     $                     crsover = .true. )
      INTEGER            NTINY
      parameter( ntiny = 11 )
      INTEGER            NL
      parameter( nl = 49 )
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d0, one = 1.0d0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, KBOT, NMIN, LLDH, LLDZ, ICTXT, NPROW, NPCOL,
     $                   myrow, mycol, hrows, hcols, ipw, nh, nb,
     $                   ii, jj, hrsrc, hcsrc, nprocs, iloc1, jloc1,
     $                   hrsrc1, hcsrc1, k, iloc2, jloc2, iloc3, jloc3,
     $                   iloc4, jloc4, hrsrc2, hcsrc2, hrsrc3, hcsrc3,
     $                   hrsrc4, hcsrc4, liwkopt
      LOGICAL            INITZ, LQUERY, WANTT, WANTZ, PAIR, BORDER
      DOUBLE PRECISION   TMP1, TMP2, TMP3, TMP4, DUM1, DUM2, DUM3,
     $                   dum4, elem1, elem4,
     $                   cs, sn, elem5, tmp, lwkopt
*     ..
*     .. Local Arrays ..
      INTEGER            DESCH2( DLEN_ )
      DOUBLE PRECISION   ELEM2( 1 ), ELEM3( 1 )
*     ..
*     .. External Functions ..
      INTEGER            PILAENVX, NUMROC, ICEIL
      LOGICAL            LSAME
      EXTERNAL           pilaenvx, lsame, numroc, iceil
*     ..
*     .. External Subroutines ..
      EXTERNAL           pdlacpy, pdlaqr1, pdlaqr0, pdlaset, pxerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dble, max, min
*     ..
*     .. Executable Statements ..
*
*     Decode and check the input parameters.
*
      info = 0
      ictxt = desch( ctxt_ )
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
      nprocs = nprow*npcol
      IF( nprow.EQ.-1 ) info = -(600+ctxt_)
      IF( info.EQ.0 ) THEN
         wantt = lsame( job, 's' )
         INITZ = LSAME( COMPZ, 'i' )
.OR.         WANTZ = INITZ  LSAME( COMPZ, 'v' )
         LLDH = DESCH( LLD_ )
         LLDZ = DESCZ( LLD_ )
         NB = DESCH( MB_ )
.EQ..OR..EQ.         LQUERY = ( LWORK-1  LIWORK-1 )
*
.NOT.         IF( LSAME( JOB, 'e.AND..NOT.' )  WANTT ) THEN
            INFO = -1
.NOT.         ELSE IF( LSAME( COMPZ, 'n.AND..NOT.' )  WANTZ ) THEN
            INFO = -2
.LT.         ELSE IF( N0 ) THEN
            INFO = -3
.LT..OR..GT.         ELSE IF( ILO1  ILOMAX( 1, N ) ) THEN
            INFO = -4
.LT..OR..GT.         ELSE IF( IHIMIN( ILO, N )  IHIN ) THEN
            INFO = -5
.NE.         ELSEIF( DESCZ( CTXT_ )DESCH( CTXT_ ) ) THEN
            INFO = -( 1000+CTXT_ )
.NE.         ELSEIF( DESCH( MB_ )DESCH( NB_ ) ) THEN
            INFO = -( 700+NB_ )
.NE.         ELSEIF( DESCZ( MB_ )DESCZ( NB_ ) ) THEN
            INFO = -( 1000+NB_ )
.NE.         ELSEIF( DESCH( MB_ )DESCZ( MB_ ) ) THEN
            INFO = -( 1000+MB_ )
.LT.         ELSEIF( DESCH( MB_ )6 ) THEN
            INFO = -( 700+NB_ )
.LT.         ELSEIF( DESCZ( MB_ )6 ) THEN
            INFO = -( 1000+MB_ )
         ELSE
            CALL CHK1MAT( N, 3, N, 3, 1, 1, DESCH, 7, INFO )
.EQ.            IF( INFO0 )
     $         CALL CHK1MAT( N, 3, N, 3, 1, 1, DESCZ, 11, INFO )
.EQ.            IF( INFO0 )
     $         CALL PCHK2MAT( N, 3, N, 3, 1, 1, DESCH, 7, N, 3, N, 3,
     $              1, 1, DESCZ, 11, 0, IWORK, IWORK, INFO )
         END IF
      END IF
*
*     Compute required workspace.
*
      CALL PDLAQR1( WANTT, WANTZ, N, ILO, IHI, H, DESCH, WR, WI,
     $     ILO, IHI, Z, DESCZ, WORK, -1, IWORK, -1, INFO )
      LWKOPT = WORK(1)
      LIWKOPT = IWORK(1)
      CALL PDLAQR0( WANTT, WANTZ, N, ILO, IHI, H, DESCH, WR, WI,
     $     ILO, IHI, Z, DESCZ, WORK, -1, IWORK, -1, INFO, 0 )
.LT.      IF( NNL ) THEN
         HROWS = NUMROC( NL, NB, MYROW, DESCH(RSRC_), NPROW )
         HCOLS = NUMROC( NL, NB, MYCOL, DESCH(CSRC_), NPCOL )
         WORK(1) = WORK(1) + DBLE(2*HROWS*HCOLS)
      END IF
      LWKOPT = MAX( LWKOPT, WORK(1) )
      LIWKOPT = MAX( LIWKOPT, IWORK(1) )
      WORK(1) = LWKOPT
      IWORK(1) = LIWKOPT
*
.NOT..AND..LT.      IF( LQUERY  LWORKINT(LWKOPT) ) THEN
         INFO = -13
.NOT..AND..LT.      ELSEIF( LQUERY  LIWORKLIWKOPT ) THEN
         INFO = -15
      END IF
*
.NE.      IF( INFO0 ) THEN
*
*        Quick return in case of invalid argument.
*
         CALL PXERBLA( ICTXT, 'pdhseqr', -INFO )
         RETURN
*
.EQ.      ELSE IF( N0 ) THEN
*
*        Quick return in case N = 0; nothing to do.
*
         RETURN
*
      ELSE IF( LQUERY ) THEN
*
*        Quick return in case of a workspace query.
*
         RETURN
*
      ELSE
*
*        Copy eigenvalues isolated by PDGEBAL.
*
         DO 10 I = 1, ILO - 1
            CALL INFOG2L( I, I, DESCH, NPROW, NPCOL, MYROW, MYCOL, II,
     $           JJ, HRSRC, HCSRC )
.EQ..AND..EQ.            IF( MYROWHRSRC  MYCOLHCSRC ) THEN
               WR( I ) = H( (JJ-1)*LLDH + II )
            ELSE
               WR( I ) = ZERO
            END IF
            WI( I ) = ZERO
   10    CONTINUE
.GT.         IF( ILO1 )
     $      CALL DGSUM2D( ICTXT, 'all', '1-tree', ILO-1, 1, WR, N, -1,
     $           -1 )
         DO 20 I = IHI + 1, N
            CALL INFOG2L( I, I, DESCH, NPROW, NPCOL, MYROW, MYCOL, II,
     $           JJ, HRSRC, HCSRC )
.EQ..AND..EQ.            IF( MYROWHRSRC  MYCOLHCSRC ) THEN
               WR( I ) = H( (JJ-1)*LLDH + II )
            ELSE
               WR( I ) = ZERO
            END IF
            WI( I ) = ZERO
   20    CONTINUE
.LT.         IF( IHIN )
     $      CALL DGSUM2D( ICTXT, 'all', '1-tree', N-IHI, 1, WR(IHI+1),
     $           N, -1, -1 )
*
*        Initialize Z, if requested.
*
         IF( INITZ )
     $      CALL PDLASET( 'a', n, n, zero, one, z, 1, 1, descz )
*
*        Quick return if possible.
*
         nprocs = nprow*npcol
         IF( ilo.EQ.ihi ) THEN
            CALL infog2l( ilo, ilo, desch, nprow, npcol, myrow,
     $           mycol, ii, jj, hrsrc, hcsrc )
            IF( myrow.EQ.hrsrc .AND. mycol.EQ.hcsrc ) THEN
               wr( ilo ) = h( (jj-1)*lldh + ii )
               IF( nprocs.GT.1 )
     $            CALL dgebs2d( ictxt, 'All', '1-Tree', 1, 1, wr(ilo),
     $                 1 )
            ELSE
               CALL dgebr2d( ictxt, 'All', '1-Tree', 1, 1, wr(ilo),
     $              1, hrsrc, hcsrc )
            END IF
            wi( ilo ) = zero
            RETURN
         END IF
*
*        PDLAQR1/PDLAQR0 crossover point.
*
         nh = ihi-ilo+1
         nmin = pilaenvx( ictxt, 12, 'PDHSEQR',
     $        job( : 1 ) // compz( : 1 ), n, ilo, ihi, lwork )
         nmin = max( ntiny, nmin )
*
*        PDLAQR0 for big matrices; PDLAQR1 for small ones.
*
         IF( (.NOT. crsover .AND. nh.GT.ntiny) .OR. nh.GT.nmin .OR.
     $        desch(rsrc_).NE.0 .OR. desch(csrc_).NE.0 ) THEN
            CALL pdlaqr0( wantt, wantz, n, ilo, ihi, h, desch, wr, wi,
     $           ilo, ihi, z, descz, work, lwork, iwork, liwork, info,
     $           0 )
            IF( info.GT.0 .AND. ( desch(rsrc_).NE.0 .OR.
     $           desch(csrc_).NE.0 ) ) THEN
*
*              A rare PDLAQR0 failure!  PDLAQR1 sometimes succeeds
*              when PDLAQR0 fails.
*
               kbot = info
               CALL pdlaqr1( wantt, wantz, n, ilo, ihi, h, desch, wr,
     $              wi, ilo, ihi, z, descz, work, lwork, iwork,
     $              liwork, info )
               info = -7777
            END IF
         ELSE
*
*           Small matrix.
*
            CALL pdlaqr1( wantt, wantz, n, ilo, ihi, h, desch, wr, wi,
     $           ilo, ihi, z, descz, work, lwork, iwork, liwork, info )
*
            IF( info.GT.0 ) THEN
*
*              A rare PDLAQR1 failure!  PDLAQR0 sometimes succeeds
*              when PDLAQR1 fails.
*
               kbot = info
*
               IF( n.GE.nl ) THEN
*
*                 Larger matrices have enough subdiagonal scratch
*                 space to call PDLAQR0 directly.
*
                  CALL pdlaqr0( wantt, wantz, n, ilo, kbot, h, desch,
     $                 wr, wi, ilo, ihi, z, descz, work, lwork,
     $                 iwork, liwork, info, 0 )
               ELSE
*
*                 Tiny matrices don't have enough subdiagonal
*                 scratch space to benefit from PDLAQR0.  Hence,
*                 tiny matrices must be copied into a larger
*                 array before calling PDLAQR0.
*
                  hrows = numroc( nl, nb, myrow, desch(rsrc_), nprow )
                  hcols = numroc( nl, nb, mycol, desch(csrc_), npcol )
                  CALL descinit( desch2, nl, nl, nb, nb, desch(rsrc_),
     $                 desch(csrc_), ictxt, max(1, hrows), info )
                  CALL pdlacpy( 'All', n, n, h, 1, 1, desch, work, 1,
     $                 1, desch2 )
                  CALL pdelset( work, n+1, n, desch2, zero )
                  CALL pdlaset( 'All', nl, nl-n, zero, zero, work, 1,
     $                 n+1, desch2 )
                  ipw = 1 + desch2(lld_)*hcols
                  CALL pdlaqr0( wantt, wantz, nl, ilo, kbot, work,
     $                 desch2, wr, wi, ilo, ihi, z, descz,
     $                 work(ipw), lwork-ipw+1, iwork,
     $                 liwork, info, 0 )
                  IF( wantt .OR. info.NE.0 )
     $               CALL pdlacpy( 'All', n, n, work, 1, 1, desch2,
     $                    h, 1, 1, desch )
               END IF
               info = -8888
            END IF
         END IF
*
*        Clear out the trash, if necessary.
*
         IF( ( wantt .OR. info.NE.0 ) .AND. n.GT.2 )
     $      CALL pdlaset( 'L', n-2, n-2, zero, zero, h, 3, 1, desch )
*
*        Force any 2-by-2 blocks to be complex conjugate pairs of
*        eigenvalues by removing false such blocks.
*
         DO 30 i = ilo, ihi-1
            CALL pdelget( 'All', ' ', tmp3, h, i+1, i, desch )
            IF( tmp3.NE.0.0d+00 ) THEN
               CALL pdelget( 'All', ' ', tmp1, h, i, i, desch )
               CALL pdelget( 'All', ' ', tmp2, h, i, i+1, desch )
               CALL pdelget( 'All', ' ', tmp4, h, i+1, i+1, desch )
               CALL dlanv2( tmp1, tmp2, tmp3, tmp4, dum1, dum2, dum3,
     $              dum4, cs, sn )
               IF( tmp3.EQ.0.0d+00 ) THEN
                  IF( wantt ) THEN
                     IF( i+2.LE.n )
     $                  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 )
                     CALL pdrot( i-1, h, 1, i, desch, 1, h, 1, i+1,
     $                    desch, 1, cs, sn, work, lwork, info )
                  END IF
                  IF( wantz ) THEN
                     CALL pdrot( n, z, 1, i, descz, 1, z, 1, i+1, descz,
     $                    1, cs, sn, work, lwork, info )
                  END IF
                  CALL pdelset( h, i, i, desch, tmp1 )
                  CALL pdelset( h, i, i+1, desch, tmp2 )
                  CALL pdelset( h, i+1, i, desch, tmp3 )
                  CALL pdelset( h, i+1, i+1, desch, tmp4 )
               END IF
            END IF
 30      CONTINUE
*
*        Read out eigenvalues: first let all the processes compute the
*        eigenvalue inside their diagonal blocks in parallel, except for
*        the eigenvalue located next to a block border. After that,
*        compute all eigenvalues located next to the block borders.
*        Finally, do a global summation over WR and WI so that all
*        processors receive the result.
*
         DO 40 k = ilo, ihi
            wr( k ) = zero
            wi( k ) = zero
 40      CONTINUE
         nb = desch( mb_ )
*
*        Loop 50: extract eigenvalues from the blocks which are not laid
*        out across a border of the processor mesh, except for those 1x1
*        blocks on the border.
*
         pair = .false.
         DO 50 k = ilo, ihi
            IF( .NOT. pair ) THEN
               border = mod( k, nb ).EQ.0 .OR. ( k.NE.1 .AND.
     $              mod( k, nb ).EQ.1 )
               IF( .NOT. border ) THEN
                  CALL infog2l( k, k, desch, nprow, npcol, myrow,
     $                 mycol, iloc1, jloc1, hrsrc1, hcsrc1 )
                  IF( myrow.EQ.hrsrc1 .AND. mycol.EQ.hcsrc1 ) THEN
                     elem1 = h((jloc1-1)*lldh+iloc1)
                     IF( k.LT.n ) THEN
                        elem3( 1 ) = h((jloc1-1)*lldh+iloc1+1)
                     ELSE
                        elem3( 1 ) = zero
                     END IF
                     IF( elem3( 1 ).NE.zero ) THEN
                        elem2( 1 ) = h((jloc1)*lldh+iloc1)
                        elem4 = h((jloc1)*lldh+iloc1+1)
                        CALL dlanv2( elem1, elem2( 1 ), elem3( 1 ),
     $                       elem4, wr( k ), wi( k ), wr( k+1 ),
     $                       wi( k+1 ), sn, cs )
                        pair = .true.
                     ELSE
                        IF( k.GT.1 ) THEN
                           tmp = h((jloc1-2)*lldh+iloc1)
                           IF( tmp.NE.zero ) THEN
                              elem1 = h((jloc1-2)*lldh+iloc1-1)
                              elem2( 1 ) = h((jloc1-1)*lldh+iloc1-1)
                              elem3( 1 ) = h((jloc1-2)*lldh+iloc1)
                              elem4 = h((jloc1-1)*lldh+iloc1)
                              CALL dlanv2( elem1, elem2( 1 ),
     $                             elem3( 1 ), elem4, wr( k-1 ),
     $                             wi( k-1 ), wr( k ), wi( k ), sn, cs )
                           ELSE
                              wr( k ) = elem1
                           END IF
                        ELSE
                           wr( k ) = elem1
                        END IF
                     END IF
                  END IF
               END IF
            ELSE
               pair = .false.
            END IF
 50      CONTINUE
*
*        Loop 60: extract eigenvalues from the blocks which are laid
*        out across a border of the processor mesh. The processors are
*        numbered as below:
*
*                        1 | 2
*                        --+--
*                        3 | 4
*
         DO 60 k = iceil(ilo,nb)*nb, ihi-1, nb
            CALL infog2l( k, k, desch, nprow, npcol, myrow, mycol,
     $           iloc1, jloc1, hrsrc1, hcsrc1 )
            CALL infog2l( k, k+1, desch, nprow, npcol, myrow, mycol,
     $           iloc2, jloc2, hrsrc2, hcsrc2 )
            CALL infog2l( k+1, k, desch, nprow, npcol, myrow, mycol,
     $           iloc3, jloc3, hrsrc3, hcsrc3 )
            CALL infog2l( k+1, k+1, desch, nprow, npcol, myrow, mycol,
     $           iloc4, jloc4, hrsrc4, hcsrc4 )
            IF( myrow.EQ.hrsrc2 .AND. mycol.EQ.hcsrc2 ) THEN
               elem2( 1 ) = h((jloc2-1)*lldh+iloc2)
               IF( hrsrc1.NE.hrsrc2 .OR. hcsrc1.NE.hcsrc2 )
     $            CALL dgesd2d( ictxt, 1, 1, elem2, 1, hrsrc1, hcsrc1)
            END IF
            IF( myrow.EQ.hrsrc3 .AND. mycol.EQ.hcsrc3 ) THEN
               elem3( 1 ) = h((jloc3-1)*lldh+iloc3)
               IF( hrsrc1.NE.hrsrc3 .OR. hcsrc1.NE.hcsrc3 )
     $            CALL dgesd2d( ictxt, 1, 1, elem3, 1, hrsrc1, hcsrc1)
            END IF
            IF( myrow.EQ.hrsrc4 .AND. mycol.EQ.hcsrc4 ) THEN
               work(1) = h((jloc4-1)*lldh+iloc4)
               IF( k+1.LT.n ) THEN
                  work(2) = h((jloc4-1)*lldh+iloc4+1)
               ELSE
                  work(2) = zero
               END IF
               IF( hrsrc1.NE.hrsrc4 .OR. hcsrc1.NE.hcsrc4 )
     $            CALL dgesd2d( ictxt, 2, 1, work, 2, hrsrc1, hcsrc1 )
            END IF
            IF( myrow.EQ.hrsrc1 .AND. mycol.EQ.hcsrc1 ) THEN
               elem1 = h((jloc1-1)*lldh+iloc1)
               IF( hrsrc1.NE.hrsrc2 .OR. hcsrc1.NE.hcsrc2 )
     $            CALL dgerv2d( ictxt, 1, 1, elem2, 1, hrsrc2, hcsrc2)
               IF( hrsrc1.NE.hrsrc3 .OR. hcsrc1.NE.hcsrc3 )
     $            CALL dgerv2d( ictxt, 1, 1, elem3, 1, hrsrc3, hcsrc3)
               IF( hrsrc1.NE.hrsrc4 .OR. hcsrc1.NE.hcsrc4 )
     $            CALL dgerv2d( ictxt, 2, 1, work, 2, hrsrc4, hcsrc4 )
               elem4 = work(1)
               elem5 = work(2)
               IF( elem5.EQ.zero ) THEN
                  IF( wr( k ).EQ.zero .AND. wi( k ).EQ.zero ) THEN
                     CALL dlanv2( elem1, elem2( 1 ), elem3( 1 ), elem4,
     $                    wr( k ), wi( k ), wr( k+1 ), wi( k+1 ), sn,
     $                    cs )
                  ELSEIF( wr( k+1 ).EQ.zero .AND. wi( k+1 ).EQ.zero )
     $                 THEN
                     wr( k+1 ) = elem4
                  END IF
               ELSEIF( wr( k ).EQ.zero .AND. wi( k ).EQ.zero )
     $              THEN
                  wr( k ) = elem1
               END IF
            END IF
 60      CONTINUE
*
         IF( nprocs.GT.1 ) THEN
            CALL dgsum2d( ictxt, 'All', ' ', ihi-ilo+1, 1, wr(ilo), n,
     $           -1, -1 )
            CALL dgsum2d( ictxt, 'All', ' ', ihi-ilo+1, 1, wi(ilo), n,
     $           -1, -1 )
         END IF
*
      END IF
*
      work(1) = lwkopt
      iwork(1) = liwkopt
      RETURN
*
*     End of PDHSEQR
*
      END