SRC/zlaein.f

   1 *> \brief \b ZLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZLAEIN + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaein.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaein.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaein.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK,
  22 *                          EPS3, SMLNUM, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       LOGICAL            NOINIT, RIGHTV
  26 *       INTEGER            INFO, LDB, LDH, N
  27 *       DOUBLE PRECISION   EPS3, SMLNUM
  28 *       COMPLEX*16         W
  29 *       ..
  30 *       .. Array Arguments ..
  31 *       DOUBLE PRECISION   RWORK( * )
  32 *       COMPLEX*16         B( LDB, * ), H( LDH, * ), V( * )
  33 *       ..
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *> ZLAEIN uses inverse iteration to find a right or left eigenvector
  42 *> corresponding to the eigenvalue W of a complex upper Hessenberg
  43 *> matrix H.
  44 *> \endverbatim
  45 *
  46 *  Arguments:
  47 *  ==========
  48 *
  49 *> \param[in] RIGHTV
  50 *> \verbatim
  51 *>          RIGHTV is LOGICAL
  52 *>          = .TRUE. : compute right eigenvector;
  53 *>          = .FALSE.: compute left eigenvector.
  54 *> \endverbatim
  55 *>
  56 *> \param[in] NOINIT
  57 *> \verbatim
  58 *>          NOINIT is LOGICAL
  59 *>          = .TRUE. : no initial vector supplied in V
  60 *>          = .FALSE.: initial vector supplied in V.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] N
  64 *> \verbatim
  65 *>          N is INTEGER
  66 *>          The order of the matrix H.  N >= 0.
  67 *> \endverbatim
  68 *>
  69 *> \param[in] H
  70 *> \verbatim
  71 *>          H is COMPLEX*16 array, dimension (LDH,N)
  72 *>          The upper Hessenberg matrix H.
  73 *> \endverbatim
  74 *>
  75 *> \param[in] LDH
  76 *> \verbatim
  77 *>          LDH is INTEGER
  78 *>          The leading dimension of the array H.  LDH >= max(1,N).
  79 *> \endverbatim
  80 *>
  81 *> \param[in] W
  82 *> \verbatim
  83 *>          W is COMPLEX*16
  84 *>          The eigenvalue of H whose corresponding right or left
  85 *>          eigenvector is to be computed.
  86 *> \endverbatim
  87 *>
  88 *> \param[in,out] V
  89 *> \verbatim
  90 *>          V is COMPLEX*16 array, dimension (N)
  91 *>          On entry, if NOINIT = .FALSE., V must contain a starting
  92 *>          vector for inverse iteration; otherwise V need not be set.
  93 *>          On exit, V contains the computed eigenvector, normalized so
  94 *>          that the component of largest magnitude has magnitude 1; here
  95 *>          the magnitude of a complex number (x,y) is taken to be
  96 *>          |x| + |y|.
  97 *> \endverbatim
  98 *>
  99 *> \param[out] B
 100 *> \verbatim
 101 *>          B is COMPLEX*16 array, dimension (LDB,N)
 102 *> \endverbatim
 103 *>
 104 *> \param[in] LDB
 105 *> \verbatim
 106 *>          LDB is INTEGER
 107 *>          The leading dimension of the array B.  LDB >= max(1,N).
 108 *> \endverbatim
 109 *>
 110 *> \param[out] RWORK
 111 *> \verbatim
 112 *>          RWORK is DOUBLE PRECISION array, dimension (N)
 113 *> \endverbatim
 114 *>
 115 *> \param[in] EPS3
 116 *> \verbatim
 117 *>          EPS3 is DOUBLE PRECISION
 118 *>          A small machine-dependent value which is used to perturb
 119 *>          close eigenvalues, and to replace zero pivots.
 120 *> \endverbatim
 121 *>
 122 *> \param[in] SMLNUM
 123 *> \verbatim
 124 *>          SMLNUM is DOUBLE PRECISION
 125 *>          A machine-dependent value close to the underflow threshold.
 126 *> \endverbatim
 127 *>
 128 *> \param[out] INFO
 129 *> \verbatim
 130 *>          INFO is INTEGER
 131 *>          = 0:  successful exit
 132 *>          = 1:  inverse iteration did not converge; V is set to the
 133 *>                last iterate.
 134 *> \endverbatim
 135 *
 136 *  Authors:
 137 *  ========
 138 *
 139 *> \author Univ. of Tennessee
 140 *> \author Univ. of California Berkeley
 141 *> \author Univ. of Colorado Denver
 142 *> \author NAG Ltd.
 143 *
 144 *> \date September 2012
 145 *
 146 *> \ingroup complex16OTHERauxiliary
 147 *
 148 *  =====================================================================
 149       SUBROUTINE ZLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK,
 150      $                   EPS3, SMLNUM, INFO )
 151 *
 152 *  -- LAPACK auxiliary routine (version 3.4.2) --
 153 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 154 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 155 *     September 2012
 156 *
 157 *     .. Scalar Arguments ..
 158       LOGICAL            NOINIT, RIGHTV
 159       INTEGER            INFO, LDB, LDH, N
 160       DOUBLE PRECISION   EPS3, SMLNUM
 161       COMPLEX*16         W
 162 *     ..
 163 *     .. Array Arguments ..
 164       DOUBLE PRECISION   RWORK( * )
 165       COMPLEX*16         B( LDB, * ), H( LDH, * ), V( * )
 166 *     ..
 167 *
 168 *  =====================================================================
 169 *
 170 *     .. Parameters ..
 171       DOUBLE PRECISION   ONE, TENTH
 172       PARAMETER          ( ONE = 1.0D+0, TENTH = 1.0D-1 )
 173       COMPLEX*16         ZERO
 174       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
 175 *     ..
 176 *     .. Local Scalars ..
 177       CHARACTER          NORMIN, TRANS
 178       INTEGER            I, IERR, ITS, J
 179       DOUBLE PRECISION   GROWTO, NRMSML, ROOTN, RTEMP, SCALE, VNORM
 180       COMPLEX*16         CDUM, EI, EJ, TEMP, X
 181 *     ..
 182 *     .. External Functions ..
 183       INTEGER            IZAMAX
 184       DOUBLE PRECISION   DZASUM, DZNRM2
 185       COMPLEX*16         ZLADIV
 186       EXTERNAL           IZAMAX, DZASUM, DZNRM2, ZLADIV
 187 *     ..
 188 *     .. External Subroutines ..
 189       EXTERNAL           ZDSCAL, ZLATRS
 190 *     ..
 191 *     .. Intrinsic Functions ..
 192       INTRINSIC          ABS, DBLE, DIMAG, MAX, SQRT
 193 *     ..
 194 *     .. Statement Functions ..
 195       DOUBLE PRECISION   CABS1
 196 *     ..
 197 *     .. Statement Function definitions ..
 198       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
 199 *     ..
 200 *     .. Executable Statements ..
 201 *
 202       INFO = 0
 203 *
 204 *     GROWTO is the threshold used in the acceptance test for an
 205 *     eigenvector.
 206 *
 207       ROOTN = SQRT( DBLE( N ) )
 208       GROWTO = TENTH / ROOTN
 209       NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
 210 *
 211 *     Form B = H - W*I (except that the subdiagonal elements are not
 212 *     stored).
 213 *
 214       DO 20 J = 1, N
 215          DO 10 I = 1, J - 1
 216             B( I, J ) = H( I, J )
 217    10    CONTINUE
 218          B( J, J ) = H( J, J ) - W
 219    20 CONTINUE
 220 *
 221       IF( NOINIT ) THEN
 222 *
 223 *        Initialize V.
 224 *
 225          DO 30 I = 1, N
 226             V( I ) = EPS3
 227    30    CONTINUE
 228       ELSE
 229 *
 230 *        Scale supplied initial vector.
 231 *
 232          VNORM = DZNRM2( N, V, 1 )
 233          CALL ZDSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), V, 1 )
 234       END IF
 235 *
 236       IF( RIGHTV ) THEN
 237 *
 238 *        LU decomposition with partial pivoting of B, replacing zero
 239 *        pivots by EPS3.
 240 *
 241          DO 60 I = 1, N - 1
 242             EI = H( I+1, I )
 243             IF( CABS1( B( I, I ) ).LT.CABS1( EI ) ) THEN
 244 *
 245 *              Interchange rows and eliminate.
 246 *
 247                X = ZLADIV( B( I, I ), EI )
 248                B( I, I ) = EI
 249                DO 40 J = I + 1, N
 250                   TEMP = B( I+1, J )
 251                   B( I+1, J ) = B( I, J ) - X*TEMP
 252                   B( I, J ) = TEMP
 253    40          CONTINUE
 254             ELSE
 255 *
 256 *              Eliminate without interchange.
 257 *
 258                IF( B( I, I ).EQ.ZERO )
 259      $            B( I, I ) = EPS3
 260                X = ZLADIV( EI, B( I, I ) )
 261                IF( X.NE.ZERO ) THEN
 262                   DO 50 J = I + 1, N
 263                      B( I+1, J ) = B( I+1, J ) - X*B( I, J )
 264    50             CONTINUE
 265                END IF
 266             END IF
 267    60    CONTINUE
 268          IF( B( N, N ).EQ.ZERO )
 269      $      B( N, N ) = EPS3
 270 *
 271          TRANS = 'N'
 272 *
 273       ELSE
 274 *
 275 *        UL decomposition with partial pivoting of B, replacing zero
 276 *        pivots by EPS3.
 277 *
 278          DO 90 J = N, 2, -1
 279             EJ = H( J, J-1 )
 280             IF( CABS1( B( J, J ) ).LT.CABS1( EJ ) ) THEN
 281 *
 282 *              Interchange columns and eliminate.
 283 *
 284                X = ZLADIV( B( J, J ), EJ )
 285                B( J, J ) = EJ
 286                DO 70 I = 1, J - 1
 287                   TEMP = B( I, J-1 )
 288                   B( I, J-1 ) = B( I, J ) - X*TEMP
 289                   B( I, J ) = TEMP
 290    70          CONTINUE
 291             ELSE
 292 *
 293 *              Eliminate without interchange.
 294 *
 295                IF( B( J, J ).EQ.ZERO )
 296      $            B( J, J ) = EPS3
 297                X = ZLADIV( EJ, B( J, J ) )
 298                IF( X.NE.ZERO ) THEN
 299                   DO 80 I = 1, J - 1
 300                      B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
 301    80             CONTINUE
 302                END IF
 303             END IF
 304    90    CONTINUE
 305          IF( B( 1, 1 ).EQ.ZERO )
 306      $      B( 1, 1 ) = EPS3
 307 *
 308          TRANS = 'C'
 309 *
 310       END IF
 311 *
 312       NORMIN = 'N'
 313       DO 110 ITS = 1, N
 314 *
 315 *        Solve U*x = scale*v for a right eigenvector
 316 *          or U**H *x = scale*v for a left eigenvector,
 317 *        overwriting x on v.
 318 *
 319          CALL ZLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, V,
 320      $                SCALE, RWORK, IERR )
 321          NORMIN = 'Y'
 322 *
 323 *        Test for sufficient growth in the norm of v.
 324 *
 325          VNORM = DZASUM( N, V, 1 )
 326          IF( VNORM.GE.GROWTO*SCALE )
 327      $      GO TO 120
 328 *
 329 *        Choose new orthogonal starting vector and try again.
 330 *
 331          RTEMP = EPS3 / ( ROOTN+ONE )
 332          V( 1 ) = EPS3
 333          DO 100 I = 2, N
 334             V( I ) = RTEMP
 335   100    CONTINUE
 336          V( N-ITS+1 ) = V( N-ITS+1 ) - EPS3*ROOTN
 337   110 CONTINUE
 338 *
 339 *     Failure to find eigenvector in N iterations.
 340 *
 341       INFO = 1
 342 *
 343   120 CONTINUE
 344 *
 345 *     Normalize eigenvector.
 346 *
 347       I = IZAMAX( N, V, 1 )
 348       CALL ZDSCAL( N, ONE / CABS1( V( I ) ), V, 1 )
 349 *
 350       RETURN
 351 *
 352 *     End of ZLAEIN
 353 *
 354       END