pshseqr.f

Go to the documentation of this file.

      SUBROUTINE pshseqr( 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( * )
      REAL               H( * ), WI( N ), WORK( * ), WR( N ), Z( * )
*     ..
*  Purpose
*  =======
*
*  PSHSEQR 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 PSGEBAL, and then passed to PSGEHRD
*          when the matrix output by PSGEBAL 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) REAL             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) REAL             array, dimension (N)
*  WI      (global output) REAL             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) REAL             array.
*          If COMPZ = 'V', on entry Z must contain the current
*          matrix Z of accumulated transformations from, e.g., PSGEHRD,
*          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) REAL             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, PSHSEQR 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: PSLAQR0 failed to converge and PSLAQR1 was called
*                   instead. This could happen. Mostly due to a bug.
*                   Please, send a bug report to the authors.
*          = -8888: PSLAQR1 failed to converge and PSLAQR0 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 PSLAQR1, PSLAQR5
*     and SLAQR6. 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 )
      REAL               ZERO, ONE
      parameter( zero = 0.0e0, one = 1.0e0 )
*     ..
*     .. 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
      REAL               TMP1, TMP2, TMP3, TMP4, DUM1, DUM2, DUM3,
     $                   dum4, elem1, elem4,
     $                   cs, sn, elem5, tmp, lwkopt
*     ..
*     .. Local Arrays ..
      INTEGER            DESCH2( DLEN_ )
      REAL               ELEM2( 1 ), ELEM3( 1 )
*     ..
*     .. External Functions ..
      INTEGER            PILAENVX, NUMROC, ICEIL
      LOGICAL            LSAME
      EXTERNAL           pilaenvx, lsame, numroc, iceil
*     ..
*     .. External Subroutines ..
      EXTERNAL           pslacpy, pslaqr1, pslaqr0, pslaset, pxerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          float, 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 pslaqr1( 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 pslaqr0( 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) + float(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, 'pshseqr', -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 PSGEBAL.
*
         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 SGSUM2D( 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 SGSUM2D( ICTXT, 'all', '1-tree', N-IHI, 1, WR(IHI+1),
     $           N, -1, -1 )
*
*        Initialize Z, if requested.
*
         IF( INITZ )
     $      CALL PSLASET( 'a', N, N, ZERO, ONE, Z, 1, 1, DESCZ )
*
*        Quick return if possible.
*
         NPROCS = NPROW*NPCOL
.EQ.         IF( ILOIHI ) THEN
            CALL INFOG2L( ILO, ILO, DESCH, NPROW, NPCOL, MYROW,
     $           MYCOL, II, JJ, HRSRC, HCSRC )
.EQ..AND..EQ.            IF( MYROWHRSRC  MYCOLHCSRC ) THEN
               WR( ILO ) = H( (JJ-1)*LLDH + II )
.GT.               IF( NPROCS1 )
     $            CALL SGEBS2D( ICTXT, 'all', '1-tree', 1, 1, WR(ILO),
     $                 1 )
            ELSE
               CALL SGEBR2D( ICTXT, 'all', '1-tree', 1, 1, WR(ILO),
     $              1, HRSRC, HCSRC )
            END IF
            WI( ILO ) = ZERO
            RETURN
         END IF
*
*        PSLAQR1/PSLAQR0 crossover point.
*
         NH = IHI-ILO+1
         NMIN = PILAENVX( ICTXT, 12, 'pshseqr',
     $        JOB( : 1 ) // COMPZ( : 1 ), N, ILO, IHI, LWORK )
         NMIN = MAX( NTINY, NMIN )
*
*        PSLAQR0 for big matrices; PSLAQR1 for small ones.
*
.NOT..AND..GT..OR..GT..OR.         IF( ( CRSOVER  NHNTINY)  NHNMIN 
.NE..OR..NE.     $        DESCH(RSRC_)0  DESCH(CSRC_)0 ) THEN
            CALL PSLAQR0( WANTT, WANTZ, N, ILO, IHI, H, DESCH, WR, WI,
     $           ILO, IHI, Z, DESCZ, WORK, LWORK, IWORK, LIWORK, INFO,
     $           0 )
.GT..AND..NE..OR.            IF( INFO0  ( DESCH(RSRC_)0 
.NE.     $           DESCH(CSRC_)0 ) ) THEN
*
*              A rare PSLAQR0 failure!  PSLAQR1 sometimes succeeds
*              when PSLAQR0 fails.
*
               KBOT = INFO
               CALL PSLAQR1( 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 PSLAQR1( WANTT, WANTZ, N, ILO, IHI, H, DESCH, WR, WI,
     $           ILO, IHI, Z, DESCZ, WORK, LWORK, IWORK, LIWORK, INFO )
*
.GT.            IF( INFO0 ) THEN
*
*              A rare PSLAQR1 failure!  PSLAQR0 sometimes succeeds
*              when PSLAQR1 fails.
*
               KBOT = INFO
*
.GE.               IF( NNL ) THEN
*
*                 Larger matrices have enough subdiagonal scratch
*                 space to call PSLAQR0 directly.
*
                  CALL PSLAQR0( 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 PSLAQR0.  Hence,
*                 tiny matrices must be copied into a larger
*                 array before calling PSLAQR0.
*
                  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 PSLACPY( 'all', N, N, H, 1, 1, DESCH, WORK, 1,
     $                 1, DESCH2 )
                  CALL PSELSET( WORK, N+1, N, DESCH2, ZERO )
                  CALL PSLASET( 'all', NL, NL-N, ZERO, ZERO, WORK, 1,
     $                 N+1, DESCH2 )
                  IPW = 1 + DESCH2(LLD_)*HCOLS
                  CALL PSLAQR0( WANTT, WANTZ, NL, ILO, KBOT, WORK,
     $                 DESCH2, WR, WI, ILO, IHI, Z, DESCZ,
     $                 WORK(IPW), LWORK-IPW+1, IWORK,
     $                 LIWORK, INFO, 0 )
.OR..NE.                  IF( WANTT  INFO0 )
     $               CALL PSLACPY( '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.
*
.OR..NE..AND..GT.         IF( ( WANTT  INFO0 )  N2 )
     $      CALL PSLASET( '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 PSELGET( 'all', ' ', tmp3, h, i+1, i, desch )
            IF( tmp3.NE.0.0e+00 ) THEN
               CALL pselget( 'All', ' ', tmp1, h, i, i, desch )
               CALL pselget( 'All', ' ', tmp2, h, i, i+1, desch )
               CALL pselget( 'All', ' ', tmp4, h, i+1, i+1, desch )
               CALL slanv2( tmp1, tmp2, tmp3, tmp4, dum1, dum2, dum3,
     $              dum4, cs, sn )
               IF( tmp3.EQ.0.0e+00 ) THEN
                  IF( wantt ) THEN
                     IF( i+2.LE.n )
     $                  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 )
                     CALL psrot( i-1, h, 1, i, desch, 1, h, 1, i+1,
     $                    desch, 1, cs, sn, work, lwork, info )
                  END IF
                  IF( wantz ) THEN
                     CALL psrot( n, z, 1, i, descz, 1, z, 1, i+1, descz,
     $                    1, cs, sn, work, lwork, info )
                  END IF
                  CALL pselset( h, i, i, desch, tmp1 )
                  CALL pselset( h, i, i+1, desch, tmp2 )
                  CALL pselset( h, i+1, i, desch, tmp3 )
                  CALL pselset( 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 slanv2( 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 slanv2( 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 sgesd2d( 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 sgesd2d( 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 sgesd2d( 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 sgerv2d( ictxt, 1, 1, elem2, 1, hrsrc2, hcsrc2)
               IF( hrsrc1.NE.hrsrc3 .OR. hcsrc1.NE.hcsrc3 )
     $            CALL sgerv2d( ictxt, 1, 1, elem3, 1, hrsrc3, hcsrc3)
               IF( hrsrc1.NE.hrsrc4 .OR. hcsrc1.NE.hcsrc4 )
     $            CALL sgerv2d( 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 slanv2( 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 sgsum2d( ictxt, 'All', ' ', ihi-ilo+1, 1, wr(ilo), n,
     $           -1, -1 )
            CALL sgsum2d( ictxt, 'All', ' ', ihi-ilo+1, 1, wi(ilo), n,
     $           -1, -1 )
         END IF
*
      END IF
*
      work(1) = lwkopt
      iwork(1) = liwkopt
      RETURN
*
*     End of PSHSEQR
*
      END