pdtrsen.f

Go to the documentation of this file.

      SUBROUTINE pdtrsen( JOB, COMPQ, SELECT, PARA, N, T, IT, JT,
     $     DESCT, Q, IQ, JQ, DESCQ, WR, WI, M, S, SEP, WORK, LWORK,
     $     IWORK, LIWORK, INFO )
*
*     Contribution from the Department of Computing Science and HPC2N,
*     Umea University, Sweden
*
*  -- ScaLAPACK computational 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 ..
      CHARACTER          COMPQ, JOB
      INTEGER            INFO, LIWORK, LWORK, M, N,
     $                   it, jt, iq, jq
      DOUBLE PRECISION   S, SEP
*     ..
*     .. Array Arguments ..
      LOGICAL            SELECT( N )
      INTEGER            PARA( 6 ), DESCT( * ), DESCQ( * ), IWORK( * )
      DOUBLE PRECISION   Q( * ), T( * ), WI( * ), WORK( * ), WR( * )
*     ..
*
*  Purpose
*  =======
*
*  PDTRSEN reorders the real Schur factorization of a real matrix
*  A = Q*T*Q**T, so that a selected cluster of eigenvalues appears
*  in the leading diagonal blocks of the upper quasi-triangular matrix
*  T, and the leading columns of Q form an orthonormal basis of the
*  corresponding right invariant subspace. The reordering is performed
*  by PDTRORD.
*
*  Optionally the routine computes the reciprocal condition numbers of
*  the cluster of eigenvalues and/or the invariant subspace. SCASY
*  library is needed for condition estimation.
*
*  T must be in Schur form (as returned by PDLAHQR), that is, block
*  upper triangular with 1-by-1 and 2-by-2 diagonal blocks.
*
*  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
*          Specifies whether condition numbers are required for the
*          cluster of eigenvalues (S) or the invariant subspace (SEP):
*          = 'N': none;
*          = 'E': for eigenvalues only (S);
*          = 'V': for invariant subspace only (SEP);
*          = 'B': for both eigenvalues and invariant subspace (S and
*                 SEP).
*
*  COMPQ   (global input) CHARACTER*1
*          = 'V': update the matrix Q of Schur vectors;
*          = 'N': do not update Q.
*
*  SELECT  (global input) LOGICAL  array, dimension (N)
*          SELECT specifies the eigenvalues in the selected cluster. To
*          select a real eigenvalue w(j), SELECT(j) must be set to
*          .TRUE.. To select a complex conjugate pair of eigenvalues
*          w(j) and w(j+1), corresponding to a 2-by-2 diagonal block,
*          either SELECT(j) or SELECT(j+1) or both must be set to
*          .TRUE.; a complex conjugate pair of eigenvalues must be
*          either both included in the cluster or both excluded.
*
*  PARA    (global input) INTEGER*6
*          Block parameters (some should be replaced by calls to
*          PILAENV and others by meaningful default values):
*          PARA(1) = maximum number of concurrent computational windows
*                    allowed in the algorithm;
*                    0 < PARA(1) <= min(NPROW,NPCOL) must hold;
*          PARA(2) = number of eigenvalues in each window;
*                    0 < PARA(2) < PARA(3) must hold;
*          PARA(3) = window size; PARA(2) < PARA(3) < DESCT(MB_)
*                    must hold;
*          PARA(4) = minimal percentage of flops required for
*                    performing matrix-matrix multiplications instead
*                    of pipelined orthogonal transformations;
*                    0 <= PARA(4) <= 100 must hold;
*          PARA(5) = width of block column slabs for row-wise
*                    application of pipelined orthogonal
*                    transformations in their factorized form;
*                    0 < PARA(5) <= DESCT(MB_) must hold.
*          PARA(6) = the maximum number of eigenvalues moved together
*                    over a process border; in practice, this will be
*                    approximately half of the cross border window size
*                    0 < PARA(6) <= PARA(2) must hold;
*
*  N       (global input) INTEGER
*          The order of the globally distributed matrix T. N >= 0.
*
*  T       (local input/output) DOUBLE PRECISION array,
*          dimension (LLD_T,LOCc(N)).
*          On entry, the local pieces of the global distributed
*          upper quasi-triangular matrix T, in Schur form. On exit, T is
*          overwritten by the local pieces of the reordered matrix T,
*          again in Schur form, with the selected eigenvalues in the
*          globally leading diagonal blocks.
*
*  IT      (global input) INTEGER
*  JT      (global input) INTEGER
*          The row and column index in the global array T indicating the
*          first column of sub( T ). IT = JT = 1 must hold.
*
*  DESCT   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the global distributed matrix T.
*
*  Q       (local input/output) DOUBLE PRECISION array,
*          dimension (LLD_Q,LOCc(N)).
*          On entry, if COMPQ = 'V', the local pieces of the global
*          distributed matrix Q of Schur vectors.
*          On exit, if COMPQ = 'V', Q has been postmultiplied by the
*          global orthogonal transformation matrix which reorders T; the
*          leading M columns of Q form an orthonormal basis for the
*          specified invariant subspace.
*          If COMPQ = 'N', Q is not referenced.
*
*  IQ      (global input) INTEGER
*  JQ      (global input) INTEGER
*          The column index in the global array Q indicating the
*          first column of sub( Q ). IQ = JQ = 1 must hold.
*
*  DESCQ   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the global distributed matrix Q.
*
*  WR      (global output) DOUBLE PRECISION array, dimension (N)
*  WI      (global output) DOUBLE PRECISION array, dimension (N)
*          The real and imaginary parts, respectively, of the reordered
*          eigenvalues of T. The eigenvalues are in principle stored in
*          the same order as on the diagonal of T, with WR(i) = T(i,i)
*          and, if T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0
*          and WI(i+1) = -WI(i).
*          Note also that if a complex eigenvalue is sufficiently
*          ill-conditioned, then its value may differ significantly
*          from its value before reordering.
*
*  M       (global output) INTEGER
*          The dimension of the specified invariant subspace.
*          0 <= M <= N.
*
*  S       (global output) DOUBLE PRECISION
*          If JOB = 'E' or 'B', S is a lower bound on the reciprocal
*          condition number for the selected cluster of eigenvalues.
*          S cannot underestimate the true reciprocal condition number
*          by more than a factor of sqrt(N). If M = 0 or N, S = 1.
*          If JOB = 'N' or 'V', S is not referenced.
*
*  SEP     (global output) DOUBLE PRECISION
*          If JOB = 'V' or 'B', SEP is the estimated reciprocal
*          condition number of the specified invariant subspace. If
*          M = 0 or N, SEP = norm(T).
*          If JOB = 'N' or 'E', SEP is not referenced.
*
*  WORK    (local workspace/output) DOUBLE PRECISION array,
*          dimension (LWORK)
*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
*  LWORK   (local input) INTEGER
*          The dimension of the array WORK.
*
*          If LWORK = -1, then a workspace query is assumed; the routine
*          only calculates the optimal size of the WORK array, returns
*          this value as the first entry of the WORK array, and no error
*          message related to LWORK is issued by PXERBLA.
*
*  IWORK   (local workspace/output) INTEGER array, dimension (LIWORK)
*
*  LIWORK  (local input) INTEGER
*          The dimension of the array IWORK.
*
*          If LIWORK = -1, then a workspace query is assumed; the
*          routine only calculates the optimal size of the IWORK array,
*          returns this value as the first entry of the IWORK array, and
*          no error message related to LIWORK is issued by PXERBLA.
*
*  INFO    (global output) INTEGER
*          = 0: successful exit
*          < 0: if INFO = -i, the i-th argument had an illegal value.
*          If the i-th argument is an array and the j-entry had
*          an illegal value, then INFO = -(i*1000+j), if the i-th
*          argument is a scalar and had an illegal value, then INFO = -i.
*          > 0: here we have several possibilites
*            *) Reordering of T failed because some eigenvalues are too
*               close to separate (the problem is very ill-conditioned);
*               T may have been partially reordered, and WR and WI
*               contain the eigenvalues in the same order as in T.
*               On exit, INFO = {the index of T where the swap failed}.
*            *) A 2-by-2 block to be reordered split into two 1-by-1
*               blocks and the second block failed to swap with an
*               adjacent block.
*               On exit, INFO = {the index of T where the swap failed}.
*            *) If INFO = N+1, there is no valid BLACS context (see the
*               BLACS documentation for details).
*            *) If INFO = N+2, the routines used in the calculation of
*               the condition numbers raised a positive warning flag
*               (see the documentation for PGESYCTD and PSYCTCON of the
*               SCASY library).
*            *) If INFO = N+3, PGESYCTD raised an input error flag;
*               please report this bug to the authors (see below).
*               If INFO = N+4, PSYCTCON raised an input error flag;
*               please report this bug to the authors (see below).
*          In a future release this subroutine may distinguish between
*          the case 1 and 2 above.
*
*  Method
*  ======
*
*  This routine performs parallel eigenvalue reordering in real Schur
*  form by parallelizing the approach proposed in [3]. The condition
*  number estimation part is performed by using techniques and code
*  from SCASY, see http://www.cs.umu.se/research/parallel/scasy.
*
*  Additional requirements
*  =======================
*
*  The following alignment requirements must hold:
*  (a) DESCT( MB_ ) = DESCT( NB_ ) = DESCQ( MB_ ) = DESCQ( NB_ )
*  (b) DESCT( RSRC_ ) = DESCQ( RSRC_ )
*  (c) DESCT( CSRC_ ) = DESCQ( CSRC_ )
*
*  All matrices must be blocked by a block factor larger than or
*  equal to two (3). This to simplify reordering across processor
*  borders in the presence of 2-by-2 blocks.
*
*  Limitations
*  ===========
*
*  This algorithm cannot work on submatrices of T and Q, i.e.,
*  IT = JT = IQ = JQ = 1 must hold. This is however no limitation
*  since PDLAHQR does not compute Schur forms of submatrices anyway.
*
*  References
*  ==========
*
*  [1] Z. Bai and J. W. Demmel; On swapping diagonal blocks in real
*      Schur form, Linear Algebra Appl., 186:73--95, 1993. Also as
*      LAPACK Working Note 54.
*
*  [2] Z. Bai, J. W. Demmel, and A. McKenney; On computing condition
*      numbers for the nonsymmetric eigenvalue problem, ACM Trans.
*      Math. Software, 19(2):202--223, 1993. Also as LAPACK Working
*      Note 13.
*
*  [3] D. Kressner; Block algorithms for reordering standard and
*      generalized Schur forms, ACM TOMS, 32(4):521-532, 2006.
*      Also LAPACK Working Note 171.
*
*  [4] R. Granat, B. Kagstrom, and D. Kressner; Parallel eigenvalue
*      reordering in real Schur form, Concurrency and Computations:
*      Practice and Experience, 21(9):1225-1250, 2009. Also as
*      LAPACK Working Note 192.
*
*  Parallel execution recommendations
*  ==================================
*
*  Use a square grid, if possible, for maximum performance. The block
*  parameters in PARA should be kept well below the data distribution
*  block size. In particular, see [3,4] for recommended settings for
*  these parameters.
*
*  In general, the parallel algorithm strives to perform as much work
*  as possible without crossing the block borders on the main block
*  diagonal.
*
*  Contributors
*  ============
*
*  Implemented by Robert Granat, Dept. of Computing Science and HPC2N,
*  Umea University, Sweden, March 2007,
*  in collaboration with Bo Kagstrom and Daniel Kressner.
*  Modified by Meiyue Shao, October 2011.
*
*  Revisions
*  =========
*
*  Please send bug-reports to granat@cs.umu.se
*
*  Keywords
*  ========
*
*  Real Schur form, eigenvalue reordering, Sylvester matrix equation
*
*  =====================================================================
*     ..
*     .. Parameters ..
      CHARACTER          TOP
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
     $                   lld_, mb_, m_, nb_, n_, rsrc_
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( TOP = '1-Tree',
     $                     block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
     $                     ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
     $                     rsrc_ = 7, csrc_ = 8, lld_ = 9,
     $                     zero = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY, WANTBH, WANTQ, WANTS, WANTSP
      INTEGER            ICOFFT12, ICTXT, IDUM1, IDUM2, IERR, ILOC1,
     $                   ipw1, iter, itt, jloc1, jtt, k, liwmin, lldt,
     $                   lldq, lwmin, myrow, mycol, n1, n2,
     $                   nb, noexsy, npcol, nprocs, nprow, space,
     $                   t12rows, t12cols, tcols, tcsrc, trows, trsrc,
     $                   wrk1, iwrk1, wrk2, iwrk2, wrk3, iwrk3
      DOUBLE PRECISION   ELEM, EST, SCALE, RNORM
*     .. Local Arrays ..
      INTEGER            DESCT12( DLEN_ ), MBNB2( 2 ), MMAX( 1 ),
     $                   MMIN( 1 )
      DOUBLE PRECISION   DPDUM1( 1 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            NUMROC
      DOUBLE PRECISION   PDLANGE
      EXTERNAL           lsame, numroc, pdlange
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, chk1mat, descinit,
     $                   igamx2d, infog2l, pdlacpy, pdtrord, pxerbla,
     $                   pchk1mat, pchk2mat
*     $                   , PGESYCTD, PSYCTCON
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max, min, sqrt
*     ..
*     .. Executable Statements ..
*
*     Get grid parameters
*
      ictxt = desct( ctxt_ )
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
      nprocs = nprow*npcol
*
*     Test if grid is O.K., i.e., the context is valid
*
      info = 0
      IF( nprow.EQ.-1 ) THEN
         info = n+1
      END IF
*
*     Check if workspace
*
      lquery = lwork.EQ.-1 .OR. liwork.EQ.-1
*
*     Test dimensions for local sanity
*
      IF( info.EQ.0 ) THEN
         CALL chk1mat( n, 5, n, 5, it, jt, desct, 9, info )
      END IF
      IF( info.EQ.0 ) THEN
         CALL chk1mat( n, 5, n, 5, iq, jq, descq, 13, info )
      END IF
*
*     Check the blocking sizes for alignment requirements
*
      IF( info.EQ.0 ) THEN
         IF( desct( mb_ ).NE.desct( nb_ ) ) info = -(1000*9 + mb_)
      END IF
      IF( info.EQ.0 ) THEN
         IF( descq( mb_ ).NE.descq( nb_ ) ) info = -(1000*13 + mb_)
      END IF
      IF( info.EQ.0 ) THEN
         IF( desct( mb_ ).NE.descq( mb_ ) ) info = -(1000*9 + mb_)
      END IF
*
*     Check the blocking sizes for minimum sizes
*
      IF( info.EQ.0 ) THEN
         IF( n.NE.desct( mb_ ) .AND. desct( mb_ ).LT.3 )
     $        info = -(1000*9 + mb_)
         IF( n.NE.descq( mb_ ) .AND. descq( mb_ ).LT.3 )
     $        info = -(1000*13 + mb_)
      END IF
*
*     Check parameters in PARA
*
      nb = desct( mb_ )
      IF( info.EQ.0 ) THEN
         IF( para(1).LT.1 .OR. para(1).GT.min(nprow,npcol) )
     $        info = -(1000 * 4 + 1)
         IF( para(2).LT.1 .OR. para(2).GE.para(3) )
     $        info = -(1000 * 4 + 2)
         IF( para(3).LT.1 .OR. para(3).GT.nb )
     $        info = -(1000 * 4 + 3)
         IF( para(4).LT.0 .OR. para(4).GT.100 )
     $        info = -(1000 * 4 + 4)
         IF( para(5).LT.1 .OR. para(5).GT.nb )
     $        info = -(1000 * 4 + 5)
         IF( para(6).LT.1 .OR. para(6).GT.para(2) )
     $        info = -(1000 * 4 + 6)
      END IF
*
*     Check requirements on IT, JT, IQ and JQ
*
      IF( info.EQ.0 ) THEN
         IF( it.NE.1 ) info = -7
         IF( jt.NE.it ) info = -8
         IF( iq.NE.1 ) info = -11
         IF( jq.NE.iq ) info = -12
      END IF
*
*     Test input parameters for global sanity
*
      IF( info.EQ.0 ) THEN
         CALL pchk1mat( n, 5, n, 5, it, jt, desct, 9, 0, idum1,
     $        idum2, info )
      END IF
      IF( info.EQ.0 ) THEN
         CALL pchk1mat( n, 5, n, 5, iq, jq, descq, 13, 0, idum1,
     $        idum2, info )
      END IF
      IF( info.EQ.0 ) THEN
         CALL pchk2mat( n, 5, n, 5, it, jt, desct, 9, n, 5, n, 5,
     $        iq, jq, descq, 13, 0, idum1, idum2, info )
      END IF
*
*     Decode and test the input parameters
*
      IF( info.EQ.0 .OR. lquery ) THEN
         wantbh = lsame( job, 'B' )
         wants = lsame( job, 'E' ) .OR. wantbh
         wantsp = lsame( job, 'v.OR.' )  WANTBH
         WANTQ = LSAME( COMPQ, 'v' )
*
.NOT.         IF( LSAME( JOB, 'n.AND..NOT..AND..NOT.' )  WANTS  WANTSP )
     $        THEN
            INFO = -1
.NOT.         ELSEIF( LSAME( COMPQ, 'n.AND..NOT.' )  WANTQ ) THEN
            INFO = -2
.LT.         ELSEIF( N0 ) THEN
            INFO = -4
         ELSE
*
*           Extract local leading dimension
*
            LLDT = DESCT( LLD_ )
            LLDQ = DESCQ( LLD_ )
*
*           Check the SELECT vector for consistency and set M to the
*           dimension of the specified invariant subspace.
*
            M = 0
            DO 10 K = 1, N
*
*              IWORK(1:N) is an integer copy of SELECT.
*
               IF( SELECT(K) ) THEN
                  IWORK(K) = 1
               ELSE
                  IWORK(K) = 0
               END IF
.LT.               IF( KN ) THEN
                  CALL INFOG2L( K+1, K, DESCT, NPROW, NPCOL,
     $                 MYROW, MYCOL, ITT, JTT, TRSRC, TCSRC )
.EQ..AND..EQ.                  IF( MYROWTRSRC  MYCOLTCSRC ) THEN
                     ELEM = T( (JTT-1)*LLDT + ITT )
.NE.                     IF( ELEMZERO ) THEN
.AND..NOT.                        IF( SELECT(K)  SELECT(K+1) ) THEN
*                           INFO = -3
                           IWORK(K+1) = 1
.NOT..AND.                        ELSEIF( SELECT(K)  SELECT(K+1) ) THEN
*                           INFO = -3
                           IWORK(K) = 1
                        END IF
                     END IF
                  END IF
               END IF
               IF( SELECT(K) ) M = M + 1
 10         CONTINUE
            MMAX( 1 ) = M
            MMIN( 1 ) = M
.GT.            IF( NPROCS1 )
     $           CALL IGAMX2D( ICTXT, 'all', TOP, 1, 1, MMAX( 1 ), 1,
     $                -1, -1, -1, -1, -1 )
.GT.            IF( NPROCS1 )
     $           CALL IGAMN2D( ICTXT, 'all', TOP, 1, 1, MMIN( 1 ), 1,
     $                -1, -1, -1, -1, -1 )
.GT.            IF( MMAX( 1 )MMIN( 1 ) ) THEN
               M = MMAX( 1 )
.GT.               IF( NPROCS1 )
     $              CALL IGAMX2D( ICTXT, 'all', TOP, N, 1, IWORK, N,
     $                   -1, -1, -1, -1, -1 )
            END IF
*
*           Set parameters for deep pipelining in parallel
*           Sylvester solver.
*
            MBNB2( 1 ) = MIN( MAX( PARA( 3 ), PARA( 2 )*2 ), NB )
            MBNB2( 2 ) = MBNB2( 1 )
*
*           Compute needed workspace
*
            N1 = M
            N2 = N - M
            IF( WANTS ) THEN
c               CALL PGESYCTD( 'Solve', 'Schur', 'Schur', 'Notranspose',
c     $              'Notranspose', -1, 'Demand', N1, N2, T, 1, 1, DESCT,
c     $              T, N1+1, N1+1, DESCT, T, 1, N1+1, DESCT, MBNB2,
c     $              WORK, -1, IWORK(N+1), -1, NOEXSY, SCALE, IERR )
               WRK1 = INT(WORK(1))
               IWRK1 = IWORK(N+1)
               WRK1 = 0
               IWRK1 = 0
            ELSE
               WRK1 = 0
               IWRK1 = 0
            END IF
*
            IF( WANTSP ) THEN
c               CALL PSYCTCON( 'Notranspose', 'Notranspose', -1,
c     $              'Demand', N1, N2, T, 1, 1, DESCT, T, N1+1, N1+1,
c     $              DESCT, MBNB2, WORK, -1, IWORK(N+1), -1, EST, ITER,
c     $              IERR )
               WRK2 = INT(WORK(1))
               IWRK2 = IWORK(N+1)
               WRK2 = 0
               IWRK2 = 0
            ELSE
               WRK2 = 0
               IWRK2 = 0
            END IF
*
            TROWS = NUMROC( N, NB, MYROW, DESCT(RSRC_), NPROW )
            TCOLS = NUMROC( N, NB, MYCOL, DESCT(CSRC_), NPCOL )
            WRK3 = N + 7*NB**2 + 2*TROWS*PARA( 3 ) + TCOLS*PARA( 3 ) +
     $           MAX( TROWS*PARA( 3 ), TCOLS*PARA( 3 ) )
            IWRK3 = 5*PARA( 1 ) + PARA(2)*PARA(3) -
     $           PARA(2) * (PARA(2) + 1 ) / 2
*
            IF( WANTSP ) THEN
               LWMIN = MAX( 1, MAX( WRK2, WRK3) )
               LIWMIN = MAX( 1, MAX( IWRK2, IWRK3 ) )+N
            ELSE IF( LSAME( JOB, 'n' ) ) THEN
               LWMIN = MAX( 1, WRK3 )
               LIWMIN = IWRK3+N
            ELSE IF( LSAME( JOB, 'e' ) ) THEN
               LWMIN = MAX( 1, MAX( WRK1, WRK3) )
               LIWMIN = MAX( 1, MAX( IWRK1, IWRK3 ) )+N
            END IF
*
.LT..AND..NOT.            IF( LWORKLWMIN  LQUERY ) THEN
               INFO = -20
.LT..AND..NOT.            ELSE IF( LIWORKLIWMIN  LQUERY ) THEN
               INFO = -22
            END IF
         END IF
      END IF
*
*     Global maximum on info
*
.GT.      IF( NPROCS1 )
     $     CALL IGAMX2D( ICTXT, 'all', TOP, 1, 1, INFO, 1, -1, -1, -1,
     $          -1, -1 )
*
*     Return if some argument is incorrect
*
.NE..AND..NOT.      IF( INFO0  LQUERY ) THEN
         M = 0
         S = ONE
         SEP = ZERO
         CALL PXERBLA( ICTXT, 'pdtrsen', -INFO )
         RETURN
      ELSEIF( LQUERY ) THEN
         WORK( 1 ) = DBLE(LWMIN)
         IWORK( 1 ) = LIWMIN
         RETURN
      END IF
*
*     Quick return if possible.
*
.EQ..OR..EQ.      IF( MN  M0 ) THEN
         IF( WANTS )
     $        S = ONE
         IF( WANTSP )
     $        SEP = PDLANGE( '1', N, N, T, IT, JT, DESCT, WORK )
         GO TO 50
      END IF
*
*     Reorder the eigenvalues.
*
      CALL PDTRORD( COMPQ, IWORK, PARA, N, T, IT, JT,
     $     DESCT, Q, IQ, JQ, DESCQ, WR, WI, M, WORK, LWORK,
     $     IWORK(N+1), LIWORK-N, INFO )
*
      IF( WANTS ) THEN
*
*        Solve Sylvester equation T11*R - R*T2 = scale*T12 for R in
*        parallel.
*
*        Copy T12 to workspace.
*
         CALL INFOG2L( 1, N1+1, DESCT, NPROW, NPCOL, MYROW,
     $        MYCOL, ILOC1, JLOC1, TRSRC, TCSRC )
         ICOFFT12 = MOD( N1, NB )
         T12ROWS = NUMROC( N1, NB, MYROW, TRSRC, NPROW )
         T12COLS = NUMROC( N2+ICOFFT12, NB, MYCOL, TCSRC, NPCOL )
         CALL DESCINIT( DESCT12, N1, N2+ICOFFT12, NB, NB, TRSRC,
     $        TCSRC, ICTXT, MAX(1,T12ROWS), IERR )
         CALL PDLACPY( 'all', N1, N2, T, 1, N1+1, DESCT, WORK,
     $        1, 1+ICOFFT12, DESCT12 )
*
*        Solve the equation to get the solution in workspace.
*
         SPACE = DESCT12( LLD_ ) * T12COLS
         IPW1 = 1 + SPACE
c         CALL PGESYCTD( 'Solve', 'Schur', 'Schur', 'Notranspose',
c     $        'Notranspose', -1, 'Demand', N1, N2, T, 1, 1, DESCT, T,
c     $        N1+1, N1+1, DESCT, WORK, 1, 1+ICOFFT12, DESCT12, MBNB2,
c     $        WORK(IPW1), LWORK-SPACE+1, IWORK(N+1), LIWORK-N, NOEXSY,
c     $        SCALE, IERR )
.LT.         IF( IERR0 ) THEN
            INFO = N+3
         ELSE
            INFO = N+2
         END IF
*
*        Estimate the reciprocal of the condition number of the cluster
*        of eigenvalues.
*
         RNORM = PDLANGE( 'frobenius', N1, N2, WORK, 1, 1+ICOFFT12,
     $        DESCT12, DPDUM1 )
.EQ.         IF( RNORMZERO ) THEN
            S = ONE
         ELSE
            S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )*
     $           SQRT( RNORM ) )
         END IF
      END IF
*
      IF( WANTSP ) THEN
*
*        Estimate sep(T11,T21) in parallel.
*
c         CALL PSYCTCON( 'Notranspose', 'Notranspose', -1, 'Demand', N1,
c     $        N2, T, 1, 1, DESCT, T, N1+1, N1+1, DESCT, MBNB2, WORK,
c     $        LWORK, IWORK(N+1), LIWORK-N, EST, ITER, IERR )
         EST = EST * SQRT(DBLE(N1*N2))
         SEP = ONE / EST
.LT.         IF( IERR0 ) THEN
            INFO = N+4
         ELSE
            INFO = N+2
         END IF
      END IF
*
*     Return to calling program.
*
 50   CONTINUE
*
      RETURN
*
*     End of PDTRSEN
*
      END
*