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' )
         wantz = initz .OR. lsame( compz, 'V' )
         lldh = desch( lld_ )
         lldz = descz( lld_ )
         nb = desch( mb_ )
         lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
*
         IF( .NOT.lsame( job, 'E' ) .AND. .NOT.wantt ) THEN
            info = -1
         ELSE IF( .NOT.lsame( compz, 'N' ) .AND. .NOT.wantz ) THEN
            info = -2
         ELSE IF( n.LT.0 ) THEN
            info = -3
         ELSE IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
            info = -4
         ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
            info = -5
         ELSEIF( descz( ctxt_ ).NE.desch( ctxt_ ) ) THEN
            info = -( 1000+ctxt_ )
         ELSEIF( desch( mb_ ).NE.desch( nb_ ) ) THEN
            info = -( 700+nb_ )
         ELSEIF( descz( mb_ ).NE.descz( nb_ ) ) THEN
            info = -( 1000+nb_ )
         ELSEIF( desch( mb_ ).NE.descz( mb_ ) ) THEN
            info = -( 1000+mb_ )
         ELSEIF( desch( mb_ ).LT.6 ) THEN
            info = -( 700+nb_ )
         ELSEIF( descz( mb_ ).LT.6 ) THEN
            info = -( 1000+mb_ )
         ELSE
            CALL chk1mat( n, 3, n, 3, 1, 1, desch, 7, info )
            IF( info.EQ.0 )
     $         CALL chk1mat( n, 3, n, 3, 1, 1, descz, 11, info )
            IF( info.EQ.0 )
     $         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 )
      IF( n.LT.nl ) 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
*
      IF( .NOT.lquery .AND. lwork.LT.int(lwkopt) ) THEN
         info = -13
      ELSEIF( .NOT.lquery .AND. liwork.LT.liwkopt ) THEN
         info = -15
      END IF
*
      IF( info.NE.0 ) THEN
*
*        Quick return in case of invalid argument.
*
         CALL pxerbla( ictxt, 'PDHSEQR', -info )
         RETURN
*
      ELSE IF( n.EQ.0 ) 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 )
            IF( myrow.EQ.hrsrc .AND. mycol.EQ.hcsrc ) THEN
               wr( i ) = h( (jj-1)*lldh + ii )
            ELSE
               wr( i ) = zero
            END IF
            wi( i ) = zero
   10    CONTINUE
         IF( ilo.GT.1 )
     $      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 )
            IF( myrow.EQ.hrsrc .AND. mycol.EQ.hcsrc ) THEN
               wr( i ) = h( (jj-1)*lldh + ii )
            ELSE
               wr( i ) = zero
            END IF
            wi( i ) = zero
   20    CONTINUE
         IF( ihi.LT.n )
     $      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 )
.EQ.               IF( TMP30.0D+00 ) THEN
                  IF( WANTT ) THEN
.LE.                     IF( I+2N )
     $                  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
.NOT.            IF(  PAIR ) THEN
.EQ..OR..NE..AND.               BORDER = MOD( K, NB )0  ( K1 
.EQ.     $              MOD( K, NB )1 )
.NOT.               IF(  BORDER ) THEN
                  CALL INFOG2L( K, K, DESCH, NPROW, NPCOL, MYROW,
     $                 MYCOL, ILOC1, JLOC1, HRSRC1, HCSRC1 )
.EQ..AND..EQ.                  IF( MYROWHRSRC1  MYCOLHCSRC1 ) THEN
                     ELEM1 = H((JLOC1-1)*LLDH+ILOC1)
.LT.                     IF( KN ) THEN
                        ELEM3( 1 ) = H((JLOC1-1)*LLDH+ILOC1+1)
                     ELSE
                        ELEM3( 1 ) = ZERO
                     END IF
.NE.                     IF( ELEM3( 1 )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
.GT.                        IF( K1 ) THEN
                           TMP = H((JLOC1-2)*LLDH+ILOC1)
.NE.                           IF( TMPZERO ) 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 )
.EQ..AND..EQ.            IF( MYROWHRSRC2  MYCOLHCSRC2 ) THEN
               ELEM2( 1 ) = H((JLOC2-1)*LLDH+ILOC2)
.NE..OR..NE.               IF( HRSRC1HRSRC2  HCSRC1HCSRC2 )
     $            CALL DGESD2D( ICTXT, 1, 1, ELEM2, 1, HRSRC1, HCSRC1)
            END IF
.EQ..AND..EQ.            IF( MYROWHRSRC3  MYCOLHCSRC3 ) THEN
               ELEM3( 1 ) = H((JLOC3-1)*LLDH+ILOC3)
.NE..OR..NE.               IF( HRSRC1HRSRC3  HCSRC1HCSRC3 )
     $            CALL DGESD2D( ICTXT, 1, 1, ELEM3, 1, HRSRC1, HCSRC1)
            END IF
.EQ..AND..EQ.            IF( MYROWHRSRC4  MYCOLHCSRC4 ) THEN
               WORK(1) = H((JLOC4-1)*LLDH+ILOC4)
.LT.               IF( K+1N ) THEN
                  WORK(2) = H((JLOC4-1)*LLDH+ILOC4+1)
               ELSE
                  WORK(2) = ZERO
               END IF
.NE..OR..NE.               IF( HRSRC1HRSRC4  HCSRC1HCSRC4 )
     $            CALL DGESD2D( ICTXT, 2, 1, WORK, 2, HRSRC1, HCSRC1 )
            END IF
.EQ..AND..EQ.            IF( MYROWHRSRC1  MYCOLHCSRC1 ) THEN
               ELEM1 = H((JLOC1-1)*LLDH+ILOC1)
.NE..OR..NE.               IF( HRSRC1HRSRC2  HCSRC1HCSRC2 )
     $            CALL DGERV2D( ICTXT, 1, 1, ELEM2, 1, HRSRC2, HCSRC2)
.NE..OR..NE.               IF( HRSRC1HRSRC3  HCSRC1HCSRC3 )
     $            CALL DGERV2D( ICTXT, 1, 1, ELEM3, 1, HRSRC3, HCSRC3)
.NE..OR..NE.               IF( HRSRC1HRSRC4  HCSRC1HCSRC4 )
     $            CALL DGERV2D( ICTXT, 2, 1, WORK, 2, HRSRC4, HCSRC4 )
               ELEM4 = WORK(1)
               ELEM5 = WORK(2)
.EQ.               IF( ELEM5ZERO ) THEN
.EQ..AND..EQ.                  IF( WR( K )ZERO  WI( K )ZERO ) THEN
                     CALL DLANV2( ELEM1, ELEM2( 1 ), ELEM3( 1 ), ELEM4,
     $                    WR( K ), WI( K ), WR( K+1 ), WI( K+1 ), SN,
     $                    CS )
.EQ..AND..EQ.                  ELSEIF( WR( K+1 )ZERO  WI( K+1 )ZERO )
     $                 THEN
                     WR( K+1 ) = ELEM4
                  END IF
.EQ..AND..EQ.               ELSEIF( WR( K )ZERO  WI( K )ZERO )
     $              THEN
                  WR( K ) = ELEM1
               END IF
            END IF
 60      CONTINUE
*
.GT.         IF( NPROCS1 ) 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