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' )
*
         IF( .NOT.lsame( job, 'N' ) .AND. .NOT.wants .AND. .NOT.wantsp )
     $        THEN
            info = -1
         ELSEIF( .NOT.lsame( compq, 'N' ) .AND. .NOT.wantq ) THEN
            info = -2
         ELSEIF( n.LT.0 ) 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
               IF( k.LT.n ) THEN
                  CALL infog2l( k+1, k, desct, nprow, npcol,
     $                 myrow, mycol, itt, jtt, trsrc, tcsrc )
                  IF( myrow.EQ.trsrc .AND. mycol.EQ.tcsrc ) THEN
                     elem = t( (jtt-1)*lldt + itt )
                     IF( elem.NE.zero ) THEN
                        IF( SELECT(k) .AND. .NOT.SELECT(k+1) ) THEN
*                           INFO = -3
                           iwork(k+1) = 1
                        ELSEIF( .NOT.SELECT(k) .AND. 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
            IF( nprocs.GT.1 )
     $           CALL igamx2d( ictxt, 'All', top, 1, 1, mmax( 1 ), 1,
     $                -1, -1, -1, -1, -1 )
            IF( nprocs.GT.1 )
     $           CALL igamn2d( ictxt, 'All', top, 1, 1, mmin( 1 ), 1,
     $                -1, -1, -1, -1, -1 )
            IF( mmax( 1 ).GT.mmin( 1 ) ) THEN
               m = mmax( 1 )
               IF( nprocs.GT.1 )
     $              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
*
            IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
               info = -20
            ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
               info = -22
            END IF
         END IF
      END IF
*
*     Global maximum on info
*
      IF( nprocs.GT.1 )
     $     CALL igamx2d( ictxt, 'All', top, 1, 1, info, 1, -1, -1, -1,
     $          -1, -1 )
*
*     Return if some argument is incorrect
*
      IF( info.NE.0 .AND. .NOT.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.
*
      IF( m.EQ.n .OR. m.EQ.0 ) 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 )
         IF( ierr.LT.0 ) 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 )
         IF( rnorm.EQ.zero ) 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
         IF( ierr.LT.0 ) THEN
            info = n+4
         ELSE
            info = n+2
         END IF
      END IF
*
*     Return to calling program.
*
 50   CONTINUE
*
      RETURN
*
*     End of PDTRSEN
*
      END
*