ScaLAPACK  2.0.2 ScaLAPACK: Scalable Linear Algebra PACKage
pztrevc.f
Go to the documentation of this file.
```00001       SUBROUTINE PZTREVC( SIDE, HOWMNY, SELECT, N, T, DESCT, VL, DESCVL,
00002      \$                    VR, DESCVR, MM, M, WORK, RWORK, INFO )
00003 *
00004 *  -- ScaLAPACK routine (version 1.7) --
00005 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
00006 *     and University of California, Berkeley.
00007 *     July 31, 2001
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          HOWMNY, SIDE
00011       INTEGER            INFO, M, MM, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       LOGICAL            SELECT( * )
00015       INTEGER            DESCT( * ), DESCVL( * ), DESCVR( * )
00016       DOUBLE PRECISION   RWORK( * )
00017       COMPLEX*16         T( * ), VL( * ), VR( * ), WORK( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  PZTREVC computes some or all of the right and/or left eigenvectors of
00024 *  a complex upper triangular matrix T in parallel.
00025 *
00026 *  The right eigenvector x and the left eigenvector y of T corresponding
00027 *  to an eigenvalue w are defined by:
00028 *
00029 *               T*x = w*x,     y'*T = w*y'
00030 *
00031 *  where y' denotes the conjugate transpose of the vector y.
00032 *
00033 *  If all eigenvectors are requested, the routine may either return the
00034 *  matrices X and/or Y of right or left eigenvectors of T, or the
00035 *  products Q*X and/or Q*Y, where Q is an input unitary
00036 *  matrix. If T was obtained from the Schur factorization of an
00037 *  original matrix A = Q*T*Q', then Q*X and Q*Y are the matrices of
00038 *  right or left eigenvectors of A.
00039 *
00040 *  Notes
00041 *  =====
00042 *
00043 *  Each global data object is described by an associated description
00044 *  vector.  This vector stores the information required to establish
00045 *  the mapping between an object element and its corresponding process
00046 *  and memory location.
00047 *
00048 *  Let A be a generic term for any 2D block cyclicly distributed array.
00049 *  Such a global array has an associated description vector DESCA.
00050 *  In the following comments, the character _ should be read as
00051 *  "of the global array".
00052 *
00053 *  NOTATION        STORED IN      EXPLANATION
00054 *  --------------- -------------- --------------------------------------
00055 *  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
00056 *                                 DTYPE_A = 1.
00057 *  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
00058 *                                 the BLACS process grid A is distribu-
00059 *                                 ted over. The context itself is glo-
00060 *                                 bal, but the handle (the integer
00061 *                                 value) may vary.
00062 *  M_A    (global) DESCA( M_ )    The number of rows in the global
00063 *                                 array A.
00064 *  N_A    (global) DESCA( N_ )    The number of columns in the global
00065 *                                 array A.
00066 *  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
00067 *                                 the rows of the array.
00068 *  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
00069 *                                 the columns of the array.
00070 *  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
00071 *                                 row of the array A is distributed.
00072 *  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
00073 *                                 first column of the array A is
00074 *                                 distributed.
00075 *  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
00076 *                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
00077 *
00078 *  Let K be the number of rows or columns of a distributed matrix,
00079 *  and assume that its process grid has dimension r x c.
00080 *  LOCr( K ) denotes the number of elements of K that a process
00081 *  would receive if K were distributed over the r processes of its
00082 *  process column.
00083 *  Similarly, LOCc( K ) denotes the number of elements of K that a
00084 *  process would receive if K were distributed over the c processes of
00085 *  its process row.
00086 *  The values of LOCr() and LOCc() may be determined via a call to the
00087 *  ScaLAPACK tool function, NUMROC:
00088 *          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
00089 *          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
00090 *  An upper bound for these quantities may be computed by:
00091 *          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
00092 *          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
00093 *
00094 *  Arguments
00095 *  =========
00096 *
00097 *  SIDE    (global input) CHARACTER*1
00098 *          = 'R':  compute right eigenvectors only;
00099 *          = 'L':  compute left eigenvectors only;
00100 *          = 'B':  compute both right and left eigenvectors.
00101 *
00102 *  HOWMNY  (global input) CHARACTER*1
00103 *          = 'A':  compute all right and/or left eigenvectors;
00104 *          = 'B':  compute all right and/or left eigenvectors,
00105 *                  and backtransform them using the input matrices
00106 *                  supplied in VR and/or VL;
00107 *          = 'S':  compute selected right and/or left eigenvectors,
00108 *                  specified by the logical array SELECT.
00109 *
00110 *  SELECT  (global input) LOGICAL array, dimension (N)
00111 *          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
00112 *          computed.
00113 *          If HOWMNY = 'A' or 'B', SELECT is not referenced.
00114 *          To select the eigenvector corresponding to the j-th
00115 *          eigenvalue, SELECT(j) must be set to .TRUE..
00116 *
00117 *  N       (global input) INTEGER
00118 *          The order of the matrix T. N >= 0.
00119 *
00120 *  T       (global input/output) COMPLEX*16 array, dimension
00121 *          (DESCT(LLD_),*)
00122 *          The upper triangular matrix T.  T is modified, but restored
00123 *          on exit.
00124 *
00125 *  DESCT   (global and local input) INTEGER array of dimension DLEN_.
00126 *          The array descriptor for the distributed matrix T.
00127 *
00128 *  VL      (global input/output) COMPLEX*16 array, dimension
00129 *          (DESCVL(LLD_),MM)
00130 *          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
00131 *          contain an N-by-N matrix Q (usually the unitary matrix Q of
00132 *          Schur vectors returned by ZHSEQR).
00133 *          On exit, if SIDE = 'L' or 'B', VL contains:
00134 *          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
00135 *          if HOWMNY = 'B', the matrix Q*Y;
00136 *          if HOWMNY = 'S', the left eigenvectors of T specified by
00137 *                           SELECT, stored consecutively in the columns
00138 *                           of VL, in the same order as their
00139 *                           eigenvalues.
00140 *          If SIDE = 'R', VL is not referenced.
00141 *
00142 *  DESCVL  (global and local input) INTEGER array of dimension DLEN_.
00143 *          The array descriptor for the distributed matrix VL.
00144 *
00145 *  VR      (global input/output) COMPLEX*16 array, dimension
00146 *          (DESCVR(LLD_),MM)
00147 *          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
00148 *          contain an N-by-N matrix Q (usually the unitary matrix Q of
00149 *          Schur vectors returned by ZHSEQR).
00150 *          On exit, if SIDE = 'R' or 'B', VR contains:
00151 *          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
00152 *          if HOWMNY = 'B', the matrix Q*X;
00153 *          if HOWMNY = 'S', the right eigenvectors of T specified by
00154 *                           SELECT, stored consecutively in the columns
00155 *                           of VR, in the same order as their
00156 *                           eigenvalues.
00157 *          If SIDE = 'L', VR is not referenced.
00158 *
00159 *  DESCVR  (global and local input) INTEGER array of dimension DLEN_.
00160 *          The array descriptor for the distributed matrix VR.
00161 *
00162 *  MM      (global input) INTEGER
00163 *          The number of columns in the arrays VL and/or VR. MM >= M.
00164 *
00165 *  M       (global output) INTEGER
00166 *          The number of columns in the arrays VL and/or VR actually
00167 *          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
00168 *          is set to N.  Each selected eigenvector occupies one
00169 *          column.
00170 *
00171 *  WORK    (local workspace) COMPLEX*16 array,
00172 *                                         dimension ( 2*DESCT(LLD_) )
00173 *          Additional workspace may be required if PZLATTRS is updated
00174 *          to use WORK.
00175 *
00176 *  RWORK   (local workspace) DOUBLE PRECISION array,
00177 *                                          dimension ( DESCT(LLD_) )
00178 *
00179 *  INFO    (global output) INTEGER
00180 *          = 0:  successful exit
00181 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00182 *
00183 *  Further Details
00184 *  ===============
00185 *
00186 *  The algorithm used in this program is basically backward (forward)
00187 *  substitution.  It is the hope that scaling would be used to make the
00188 *  the code robust against possible overflow.  But scaling has not yet
00189 *  been implemented in PZLATTRS which is called by this routine to solve
00190 *  the triangular systems.  PZLATTRS just calls PZTRSV.
00191 *
00192 *  Each eigenvector is normalized so that the element of largest
00193 *  magnitude has magnitude 1; here the magnitude of a complex number
00194 *  (x,y) is taken to be |x| + |y|.
00195 *
00196 *  Further Details
00197 *  ===============
00198 *
00199 *  Implemented by Mark R. Fahey, June, 2000
00200 *
00201 *  =====================================================================
00202 *
00203 *     .. Parameters ..
00204       DOUBLE PRECISION   ZERO, ONE
00205       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00206       COMPLEX*16         CZERO, CONE
00207       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00208      \$                   CONE = ( 1.0D+0, 0.0D+0 ) )
00209       INTEGER            BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
00210      \$                   MB_, NB_, RSRC_, CSRC_, LLD_
00211       PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
00212      \$                   CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
00213      \$                   RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
00214 *     ..
00215 *     .. Local Scalars ..
00216       LOGICAL            ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV
00217       INTEGER            CONTXT, CSRC, I, ICOL, II, IROW, IS, ITMP1,
00218      \$                   ITMP2, J, K, KI, LDT, LDVL, LDVR, LDW, MB,
00219      \$                   MYCOL, MYROW, NB, NPCOL, NPROW, RSRC
00220       REAL               SELF
00221       DOUBLE PRECISION   OVFL, REMAXD, SCALE, SMIN, SMLNUM, ULP, UNFL
00222       COMPLEX*16         CDUM, REMAXC, SHIFT
00223 *     ..
00224 *     .. Local Arrays ..
00225       INTEGER            DESCW( DLEN_ )
00226 *     ..
00227 *     .. External Functions ..
00228       LOGICAL            LSAME
00229       DOUBLE PRECISION   PDLAMCH
00230       EXTERNAL           LSAME, PDLAMCH
00231 *     ..
00232 *     .. External Subroutines ..
00233       EXTERNAL           BLACS_GRIDINFO, DESCINIT, DGSUM2D, IGAMN2D,
00234      \$                   INFOG2L, PDLABAD, PDZASUM, PXERBLA, PZAMAX,
00235      \$                   PZCOPY, PZDSCAL, PZGEMV, PZLASET, PZLATTRS,
00236      \$                   ZGSUM2D
00237 *     ..
00238 *     .. Intrinsic Functions ..
00239       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX
00240 *     ..
00241 *     .. Statement Functions ..
00242       DOUBLE PRECISION   CABS1
00243 *     ..
00244 *     .. Statement Function definitions ..
00245       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
00246 *     ..
00247 *     .. Executable Statements ..
00248 *
00249 *       This is just to keep ftnchek happy
00250       IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_*
00251      \$    RSRC_.LT.0 )RETURN
00252 *
00253       CONTXT = DESCT( CTXT_ )
00254       RSRC = DESCT( RSRC_ )
00255       CSRC = DESCT( CSRC_ )
00256       MB = DESCT( MB_ )
00257       NB = DESCT( NB_ )
00258       LDT = DESCT( LLD_ )
00259       LDW = LDT
00260       LDVR = DESCVR( LLD_ )
00261       LDVL = DESCVL( LLD_ )
00262 *
00263       CALL BLACS_GRIDINFO( CONTXT, NPROW, NPCOL, MYROW, MYCOL )
00264       SELF = MYROW*NPCOL + MYCOL
00265 *
00266 *     Decode and test the input parameters
00267 *
00268       BOTHV = LSAME( SIDE, 'B' )
00269       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
00270       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
00271 *
00272       ALLV = LSAME( HOWMNY, 'A' )
00273       OVER = LSAME( HOWMNY, 'B' ) .OR. LSAME( HOWMNY, 'O' )
00274       SOMEV = LSAME( HOWMNY, 'S' )
00275 *
00276 *     Set M to the number of columns required to store the selected
00277 *     eigenvectors.
00278 *
00279       IF( SOMEV ) THEN
00280          M = 0
00281          DO 10 J = 1, N
00282             IF( SELECT( J ) )
00283      \$         M = M + 1
00284    10    CONTINUE
00285       ELSE
00286          M = N
00287       END IF
00288 *
00289       INFO = 0
00290       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
00291          INFO = -1
00292       ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
00293          INFO = -2
00294       ELSE IF( N.LT.0 ) THEN
00295          INFO = -4
00296       ELSE IF( MM.LT.M ) THEN
00297          INFO = -11
00298       END IF
00299       CALL IGAMN2D( CONTXT, 'ALL', ' ', 1, 1, INFO, 1, ITMP1, ITMP2, -1,
00300      \$              -1, -1 )
00301       IF( INFO.LT.0 ) THEN
00302          CALL PXERBLA( CONTXT, 'PZTREVC', -INFO )
00303          RETURN
00304       END IF
00305 *
00306 *     Quick return if possible.
00307 *
00308       IF( N.EQ.0 )
00309      \$   RETURN
00310 *
00311 *     Set the constants to control overflow.
00312 *
00313       UNFL = PDLAMCH( CONTXT, 'Safe minimum' )
00314       OVFL = ONE / UNFL
00315       CALL PDLABAD( CONTXT, UNFL, OVFL )
00316       ULP = PDLAMCH( CONTXT, 'Precision' )
00317       SMLNUM = UNFL*( N / ULP )
00318 *
00319 *     Store the diagonal elements of T in working array WORK( LDW+1 ).
00320 *
00321       DO 20 I = 1, N
00322          CALL INFOG2L( I, I, DESCT, NPROW, NPCOL, MYROW, MYCOL, IROW,
00323      \$                 ICOL, ITMP1, ITMP2 )
00324          IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN
00325             WORK( LDW+IROW ) = T( ( ICOL-1 )*LDT+IROW )
00326          END IF
00327    20 CONTINUE
00328 *
00329 *     Compute 1-norm of each column of strictly upper triangular
00330 *     part of T to control overflow in triangular solver.  Computed,
00331 *     but not used.  For use in PZLATTRS.
00332 *
00333       RWORK( 1 ) = ZERO
00334       DO 30 J = 2, N
00335          CALL PDZASUM( J-1, RWORK( J ), T, 1, J, DESCT, 1 )
00336    30 CONTINUE
00337 *     I replicate the norms in RWORK.  Should they be distributed
00338 *     over the process rows?
00339       CALL DGSUM2D( CONTXT, 'Row', ' ', N, 1, RWORK, N, -1, -1 )
00340 *
00341       IF( RIGHTV ) THEN
00342 *
00343 *        Compute right eigenvectors.
00344 *
00345 *        Need to set the distribution pattern of WORK
00346 *
00347          CALL DESCINIT( DESCW, N, 1, NB, 1, RSRC, CSRC, CONTXT, LDW,
00348      \$                  INFO )
00349 *
00350          IS = M
00351          DO 70 KI = N, 1, -1
00352 *
00353             IF( SOMEV ) THEN
00354                IF( .NOT.SELECT( KI ) )
00355      \$            GO TO 70
00356             END IF
00357 *
00358             SMIN = ZERO
00359             SHIFT = CZERO
00360             CALL INFOG2L( KI, KI, DESCT, NPROW, NPCOL, MYROW, MYCOL,
00361      \$                    IROW, ICOL, ITMP1, ITMP2 )
00362             IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN
00363                SHIFT = T( ( ICOL-1 )*LDT+IROW )
00364                SMIN = MAX( ULP*( CABS1( SHIFT ) ), SMLNUM )
00365             END IF
00366             CALL DGSUM2D( CONTXT, 'ALL', ' ', 1, 1, SMIN, 1, -1, -1 )
00367             CALL ZGSUM2D( CONTXT, 'ALL', ' ', 1, 1, SHIFT, 1, -1, -1 )
00368 *
00369             CALL INFOG2L( 1, 1, DESCW, NPROW, NPCOL, MYROW, MYCOL, IROW,
00370      \$                    ICOL, ITMP1, ITMP2 )
00371             IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN
00372                WORK( 1 ) = CONE
00373             END IF
00374 *
00375 *           Form right-hand side.  Distribute rhs onto first column
00376 *           of processor grid.
00377 *
00378             IF( KI.GT.1 ) THEN
00379                CALL PZCOPY( KI-1, T, 1, KI, DESCT, 1, WORK, 1, 1, DESCW,
00380      \$                      1 )
00381             END IF
00382             DO 40 K = 1, KI - 1
00383                CALL INFOG2L( K, 1, DESCW, NPROW, NPCOL, MYROW, MYCOL,
00384      \$                       IROW, ICOL, ITMP1, ITMP2 )
00385                IF( MYROW.EQ.ITMP1 .AND. MYCOL.EQ.ITMP2 ) THEN
00386                   WORK( IROW ) = -WORK( IROW )
00387                END IF
00388    40       CONTINUE
00389 *
00390 *           Solve the triangular system:
00391 *              (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK.
00392 *
00393             DO 50 K = 1, KI - 1
00394                CALL INFOG2L( K, K, DESCT, NPROW, NPCOL, MYROW, MYCOL,
00395      \$                       IROW, ICOL, ITMP1, ITMP2 )
00396                IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN
00397                   T( ( ICOL-1 )*LDT+IROW ) = T( ( ICOL-1 )*LDT+IROW ) -
00398      \$               SHIFT
00399                   IF( CABS1( T( ( ICOL-1 )*LDT+IROW ) ).LT.SMIN ) THEN
00400                      T( ( ICOL-1 )*LDT+IROW ) = DCMPLX( SMIN )
00401                   END IF
00402                END IF
00403    50       CONTINUE
00404 *
00405             IF( KI.GT.1 ) THEN
00406                CALL PZLATTRS( 'Upper', 'No transpose', 'Non-unit', 'Y',
00407      \$                        KI-1, T, 1, 1, DESCT, WORK, 1, 1, DESCW,
00408      \$                        SCALE, RWORK, INFO )
00409                CALL INFOG2L( KI, 1, DESCW, NPROW, NPCOL, MYROW, MYCOL,
00410      \$                       IROW, ICOL, ITMP1, ITMP2 )
00411                IF( MYROW.EQ.ITMP1 .AND. MYCOL.EQ.ITMP2 ) THEN
00412                   WORK( IROW ) = DCMPLX( SCALE )
00413                END IF
00414             END IF
00415 *
00416 *           Copy the vector x or Q*x to VR and normalize.
00417 *
00418             IF( .NOT.OVER ) THEN
00419                CALL PZCOPY( KI, WORK, 1, 1, DESCW, 1, VR, 1, IS, DESCVR,
00420      \$                      1 )
00421 *
00422                CALL PZAMAX( KI, REMAXC, II, VR, 1, IS, DESCVR, 1 )
00423                REMAXD = ONE / MAX( CABS1( REMAXC ), UNFL )
00424                CALL PZDSCAL( KI, REMAXD, VR, 1, IS, DESCVR, 1 )
00425 *
00426                CALL PZLASET( ' ', N-KI, 1, CZERO, CZERO, VR, KI+1, IS,
00427      \$                       DESCVR )
00428             ELSE
00429                IF( KI.GT.1 )
00430      \$            CALL PZGEMV( 'N', N, KI-1, CONE, VR, 1, 1, DESCVR,
00431      \$                         WORK, 1, 1, DESCW, 1, DCMPLX( SCALE ),
00432      \$                         VR, 1, KI, DESCVR, 1 )
00433 *
00434                CALL PZAMAX( N, REMAXC, II, VR, 1, KI, DESCVR, 1 )
00435                REMAXD = ONE / MAX( CABS1( REMAXC ), UNFL )
00436                CALL PZDSCAL( N, REMAXD, VR, 1, KI, DESCVR, 1 )
00437             END IF
00438 *
00439 *           Set back the original diagonal elements of T.
00440 *
00441             DO 60 K = 1, KI - 1
00442                CALL INFOG2L( K, K, DESCT, NPROW, NPCOL, MYROW, MYCOL,
00443      \$                       IROW, ICOL, ITMP1, ITMP2 )
00444                IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN
00445                   T( ( ICOL-1 )*LDT+IROW ) = WORK( LDW+IROW )
00446                END IF
00447    60       CONTINUE
00448 *
00449             IS = IS - 1
00450    70    CONTINUE
00451       END IF
00452 *
00453       IF( LEFTV ) THEN
00454 *
00455 *        Compute left eigenvectors.
00456 *
00457 *        Need to set the distribution pattern of WORK
00458 *
00459          CALL DESCINIT( DESCW, N, 1, MB, 1, RSRC, CSRC, CONTXT, LDW,
00460      \$                  INFO )
00461 *
00462          IS = 1
00463          DO 110 KI = 1, N
00464 *
00465             IF( SOMEV ) THEN
00466                IF( .NOT.SELECT( KI ) )
00467      \$            GO TO 110
00468             END IF
00469 *
00470             SMIN = ZERO
00471             SHIFT = CZERO
00472             CALL INFOG2L( KI, KI, DESCT, NPROW, NPCOL, MYROW, MYCOL,
00473      \$                    IROW, ICOL, ITMP1, ITMP2 )
00474             IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN
00475                SHIFT = T( ( ICOL-1 )*LDT+IROW )
00476                SMIN = MAX( ULP*( CABS1( SHIFT ) ), SMLNUM )
00477             END IF
00478             CALL DGSUM2D( CONTXT, 'ALL', ' ', 1, 1, SMIN, 1, -1, -1 )
00479             CALL ZGSUM2D( CONTXT, 'ALL', ' ', 1, 1, SHIFT, 1, -1, -1 )
00480 *
00481             CALL INFOG2L( N, 1, DESCW, NPROW, NPCOL, MYROW, MYCOL, IROW,
00482      \$                    ICOL, ITMP1, ITMP2 )
00483             IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN
00484                WORK( IROW ) = CONE
00485             END IF
00486 *
00487 *           Form right-hand side.
00488 *
00489             IF( KI.LT.N ) THEN
00490                CALL PZCOPY( N-KI, T, KI, KI+1, DESCT, N, WORK, KI+1, 1,
00491      \$                      DESCW, 1 )
00492             END IF
00493             DO 80 K = KI + 1, N
00494                CALL INFOG2L( K, 1, DESCW, NPROW, NPCOL, MYROW, MYCOL,
00495      \$                       IROW, ICOL, ITMP1, ITMP2 )
00496                IF( MYROW.EQ.ITMP1 .AND. MYCOL.EQ.ITMP2 ) THEN
00497                   WORK( IROW ) = -DCONJG( WORK( IROW ) )
00498                END IF
00499    80       CONTINUE
00500 *
00501 *           Solve the triangular system:
00502 *              (T(KI+1:N,KI+1:N) - T(KI,KI))'*X = SCALE*WORK.
00503 *
00504             DO 90 K = KI + 1, N
00505                CALL INFOG2L( K, K, DESCT, NPROW, NPCOL, MYROW, MYCOL,
00506      \$                       IROW, ICOL, ITMP1, ITMP2 )
00507                IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN
00508                   T( ( ICOL-1 )*LDT+IROW ) = T( ( ICOL-1 )*LDT+IROW ) -
00509      \$               SHIFT
00510                   IF( CABS1( T( ( ICOL-1 )*LDT+IROW ) ).LT.SMIN )
00511      \$               T( ( ICOL-1 )*LDT+IROW ) = DCMPLX( SMIN )
00512                END IF
00513    90       CONTINUE
00514 *
00515             IF( KI.LT.N ) THEN
00516                CALL PZLATTRS( 'Upper', 'Conjugate transpose', 'Nonunit',
00517      \$                        'Y', N-KI, T, KI+1, KI+1, DESCT, WORK,
00518      \$                        KI+1, 1, DESCW, SCALE, RWORK, INFO )
00519                CALL INFOG2L( KI, 1, DESCW, NPROW, NPCOL, MYROW, MYCOL,
00520      \$                       IROW, ICOL, ITMP1, ITMP2 )
00521                IF( MYROW.EQ.ITMP1 .AND. MYCOL.EQ.ITMP2 ) THEN
00522                   WORK( IROW ) = DCMPLX( SCALE )
00523                END IF
00524             END IF
00525 *
00526 *           Copy the vector x or Q*x to VL and normalize.
00527 *
00528             IF( .NOT.OVER ) THEN
00529                CALL PZCOPY( N-KI+1, WORK, KI, 1, DESCW, 1, VL, KI, IS,
00530      \$                      DESCVL, 1 )
00531 *
00532                CALL PZAMAX( N-KI+1, REMAXC, II, VL, KI, IS, DESCVL, 1 )
00533                REMAXD = ONE / MAX( CABS1( REMAXC ), UNFL )
00534                CALL PZDSCAL( N-KI+1, REMAXD, VL, KI, IS, DESCVL, 1 )
00535 *
00536                CALL PZLASET( ' ', KI-1, 1, CZERO, CZERO, VL, 1, IS,
00537      \$                       DESCVL )
00538             ELSE
00539                IF( KI.LT.N )
00540      \$            CALL PZGEMV( 'N', N, N-KI, CONE, VL, 1, KI+1, DESCVL,
00541      \$                         WORK, KI+1, 1, DESCW, 1, DCMPLX( SCALE ),
00542      \$                         VL, 1, KI, DESCVL, 1 )
00543 *
00544                CALL PZAMAX( N, REMAXC, II, VL, 1, KI, DESCVL, 1 )
00545                REMAXD = ONE / MAX( CABS1( REMAXC ), UNFL )
00546                CALL PZDSCAL( N, REMAXD, VL, 1, KI, DESCVL, 1 )
00547             END IF
00548 *
00549 *           Set back the original diagonal elements of T.
00550 *
00551             DO 100 K = KI + 1, N
00552                CALL INFOG2L( K, K, DESCT, NPROW, NPCOL, MYROW, MYCOL,
00553      \$                       IROW, ICOL, ITMP1, ITMP2 )
00554                IF( ( MYROW.EQ.ITMP1 ) .AND. ( MYCOL.EQ.ITMP2 ) ) THEN
00555                   T( ( ICOL-1 )*LDT+IROW ) = WORK( LDW+IROW )
00556                END IF
00557   100       CONTINUE
00558 *
00559             IS = IS + 1
00560   110    CONTINUE
00561       END IF
00562 *
00563       RETURN
00564 *
00565 *     End of PZTREVC
00566 *
00567       END
```