SRC/zstein.f

   1 *> \brief \b ZSTEIN
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZSTEIN + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstein.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstein.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstein.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
  22 *                          IWORK, IFAIL, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            INFO, LDZ, M, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       INTEGER            IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
  29 *      $                   IWORK( * )
  30 *       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
  31 *       COMPLEX*16         Z( LDZ, * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> ZSTEIN computes the eigenvectors of a real symmetric tridiagonal
  41 *> matrix T corresponding to specified eigenvalues, using inverse
  42 *> iteration.
  43 *>
  44 *> The maximum number of iterations allowed for each eigenvector is
  45 *> specified by an internal parameter MAXITS (currently set to 5).
  46 *>
  47 *> Although the eigenvectors are real, they are stored in a complex
  48 *> array, which may be passed to ZUNMTR or ZUPMTR for back
  49 *> transformation to the eigenvectors of a complex Hermitian matrix
  50 *> which was reduced to tridiagonal form.
  51 *>
  52 *> \endverbatim
  53 *
  54 *  Arguments:
  55 *  ==========
  56 *
  57 *> \param[in] N
  58 *> \verbatim
  59 *>          N is INTEGER
  60 *>          The order of the matrix.  N >= 0.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] D
  64 *> \verbatim
  65 *>          D is DOUBLE PRECISION array, dimension (N)
  66 *>          The n diagonal elements of the tridiagonal matrix T.
  67 *> \endverbatim
  68 *>
  69 *> \param[in] E
  70 *> \verbatim
  71 *>          E is DOUBLE PRECISION array, dimension (N-1)
  72 *>          The (n-1) subdiagonal elements of the tridiagonal matrix
  73 *>          T, stored in elements 1 to N-1.
  74 *> \endverbatim
  75 *>
  76 *> \param[in] M
  77 *> \verbatim
  78 *>          M is INTEGER
  79 *>          The number of eigenvectors to be found.  0 <= M <= N.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] W
  83 *> \verbatim
  84 *>          W is DOUBLE PRECISION array, dimension (N)
  85 *>          The first M elements of W contain the eigenvalues for
  86 *>          which eigenvectors are to be computed.  The eigenvalues
  87 *>          should be grouped by split-off block and ordered from
  88 *>          smallest to largest within the block.  ( The output array
  89 *>          W from DSTEBZ with ORDER = 'B' is expected here. )
  90 *> \endverbatim
  91 *>
  92 *> \param[in] IBLOCK
  93 *> \verbatim
  94 *>          IBLOCK is INTEGER array, dimension (N)
  95 *>          The submatrix indices associated with the corresponding
  96 *>          eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
  97 *>          the first submatrix from the top, =2 if W(i) belongs to
  98 *>          the second submatrix, etc.  ( The output array IBLOCK
  99 *>          from DSTEBZ is expected here. )
 100 *> \endverbatim
 101 *>
 102 *> \param[in] ISPLIT
 103 *> \verbatim
 104 *>          ISPLIT is INTEGER array, dimension (N)
 105 *>          The splitting points, at which T breaks up into submatrices.
 106 *>          The first submatrix consists of rows/columns 1 to
 107 *>          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
 108 *>          through ISPLIT( 2 ), etc.
 109 *>          ( The output array ISPLIT from DSTEBZ is expected here. )
 110 *> \endverbatim
 111 *>
 112 *> \param[out] Z
 113 *> \verbatim
 114 *>          Z is COMPLEX*16 array, dimension (LDZ, M)
 115 *>          The computed eigenvectors.  The eigenvector associated
 116 *>          with the eigenvalue W(i) is stored in the i-th column of
 117 *>          Z.  Any vector which fails to converge is set to its current
 118 *>          iterate after MAXITS iterations.
 119 *>          The imaginary parts of the eigenvectors are set to zero.
 120 *> \endverbatim
 121 *>
 122 *> \param[in] LDZ
 123 *> \verbatim
 124 *>          LDZ is INTEGER
 125 *>          The leading dimension of the array Z.  LDZ >= max(1,N).
 126 *> \endverbatim
 127 *>
 128 *> \param[out] WORK
 129 *> \verbatim
 130 *>          WORK is DOUBLE PRECISION array, dimension (5*N)
 131 *> \endverbatim
 132 *>
 133 *> \param[out] IWORK
 134 *> \verbatim
 135 *>          IWORK is INTEGER array, dimension (N)
 136 *> \endverbatim
 137 *>
 138 *> \param[out] IFAIL
 139 *> \verbatim
 140 *>          IFAIL is INTEGER array, dimension (M)
 141 *>          On normal exit, all elements of IFAIL are zero.
 142 *>          If one or more eigenvectors fail to converge after
 143 *>          MAXITS iterations, then their indices are stored in
 144 *>          array IFAIL.
 145 *> \endverbatim
 146 *>
 147 *> \param[out] INFO
 148 *> \verbatim
 149 *>          INFO is INTEGER
 150 *>          = 0: successful exit
 151 *>          < 0: if INFO = -i, the i-th argument had an illegal value
 152 *>          > 0: if INFO = i, then i eigenvectors failed to converge
 153 *>               in MAXITS iterations.  Their indices are stored in
 154 *>               array IFAIL.
 155 *> \endverbatim
 156 *
 157 *> \par Internal Parameters:
 158 *  =========================
 159 *>
 160 *> \verbatim
 161 *>  MAXITS  INTEGER, default = 5
 162 *>          The maximum number of iterations performed.
 163 *>
 164 *>  EXTRA   INTEGER, default = 2
 165 *>          The number of iterations performed after norm growth
 166 *>          criterion is satisfied, should be at least 1.
 167 *> \endverbatim
 168 *
 169 *  Authors:
 170 *  ========
 171 *
 172 *> \author Univ. of Tennessee
 173 *> \author Univ. of California Berkeley
 174 *> \author Univ. of Colorado Denver
 175 *> \author NAG Ltd.
 176 *
 177 *> \date November 2015
 178 *
 179 *> \ingroup complex16OTHERcomputational
 180 *
 181 *  =====================================================================
 182       SUBROUTINE ZSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
 183      $                   IWORK, IFAIL, INFO )
 184 *
 185 *  -- LAPACK computational routine (version 3.6.0) --
 186 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 187 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 188 *     November 2015
 189 *
 190 *     .. Scalar Arguments ..
 191       INTEGER            INFO, LDZ, M, N
 192 *     ..
 193 *     .. Array Arguments ..
 194       INTEGER            IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
 195      $                   IWORK( * )
 196       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
 197       COMPLEX*16         Z( LDZ, * )
 198 *     ..
 199 *
 200 * =====================================================================
 201 *
 202 *     .. Parameters ..
 203       COMPLEX*16         CZERO, CONE
 204       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
 205      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
 206       DOUBLE PRECISION   ZERO, ONE, TEN, ODM3, ODM1
 207       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1,
 208      $                   ODM3 = 1.0D-3, ODM1 = 1.0D-1 )
 209       INTEGER            MAXITS, EXTRA
 210       PARAMETER          ( MAXITS = 5, EXTRA = 2 )
 211 *     ..
 212 *     .. Local Scalars ..
 213       INTEGER            B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
 214      $                   INDRV2, INDRV3, INDRV4, INDRV5, ITS, J, J1,
 215      $                   JBLK, JMAX, JR, NBLK, NRMCHK
 216       DOUBLE PRECISION   DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
 217      $                   SCL, SEP, TOL, XJ, XJM, ZTR
 218 *     ..
 219 *     .. Local Arrays ..
 220       INTEGER            ISEED( 4 )
 221 *     ..
 222 *     .. External Functions ..
 223       INTEGER            IDAMAX
 224       DOUBLE PRECISION   DLAMCH, DNRM2
 225       EXTERNAL           IDAMAX, DLAMCH, DNRM2
 226 *     ..
 227 *     .. External Subroutines ..
 228       EXTERNAL           DCOPY, DLAGTF, DLAGTS, DLARNV, DSCAL, XERBLA
 229 *     ..
 230 *     .. Intrinsic Functions ..
 231       INTRINSIC          ABS, DBLE, DCMPLX, MAX, SQRT
 232 *     ..
 233 *     .. Executable Statements ..
 234 *
 235 *     Test the input parameters.
 236 *
 237       INFO = 0
 238       DO 10 I = 1, M
 239          IFAIL( I ) = 0
 240    10 CONTINUE
 241 *
 242       IF( N.LT.0 ) THEN
 243          INFO = -1
 244       ELSE IF( M.LT.0 .OR. M.GT.N ) THEN
 245          INFO = -4
 246       ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN
 247          INFO = -9
 248       ELSE
 249          DO 20 J = 2, M
 250             IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN
 251                INFO = -6
 252                GO TO 30
 253             END IF
 254             IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) )
 255      $           THEN
 256                INFO = -5
 257                GO TO 30
 258             END IF
 259    20    CONTINUE
 260    30    CONTINUE
 261       END IF
 262 *
 263       IF( INFO.NE.0 ) THEN
 264          CALL XERBLA( 'ZSTEIN', -INFO )
 265          RETURN
 266       END IF
 267 *
 268 *     Quick return if possible
 269 *
 270       IF( N.EQ.0 .OR. M.EQ.0 ) THEN
 271          RETURN
 272       ELSE IF( N.EQ.1 ) THEN
 273          Z( 1, 1 ) = CONE
 274          RETURN
 275       END IF
 276 *
 277 *     Get machine constants.
 278 *
 279       EPS = DLAMCH( 'Precision' )
 280 *
 281 *     Initialize seed for random number generator DLARNV.
 282 *
 283       DO 40 I = 1, 4
 284          ISEED( I ) = 1
 285    40 CONTINUE
 286 *
 287 *     Initialize pointers.
 288 *
 289       INDRV1 = 0
 290       INDRV2 = INDRV1 + N
 291       INDRV3 = INDRV2 + N
 292       INDRV4 = INDRV3 + N
 293       INDRV5 = INDRV4 + N
 294 *
 295 *     Compute eigenvectors of matrix blocks.
 296 *
 297       J1 = 1
 298       DO 180 NBLK = 1, IBLOCK( M )
 299 *
 300 *        Find starting and ending indices of block nblk.
 301 *
 302          IF( NBLK.EQ.1 ) THEN
 303             B1 = 1
 304          ELSE
 305             B1 = ISPLIT( NBLK-1 ) + 1
 306          END IF
 307          BN = ISPLIT( NBLK )
 308          BLKSIZ = BN - B1 + 1
 309          IF( BLKSIZ.EQ.1 )
 310      $      GO TO 60
 311          GPIND = J1
 312 *
 313 *        Compute reorthogonalization criterion and stopping criterion.
 314 *
 315          ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) )
 316          ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) )
 317          DO 50 I = B1 + 1, BN - 1
 318             ONENRM = MAX( ONENRM, ABS( D( I ) )+ABS( E( I-1 ) )+
 319      $               ABS( E( I ) ) )
 320    50    CONTINUE
 321          ORTOL = ODM3*ONENRM
 322 *
 323          DTPCRT = SQRT( ODM1 / BLKSIZ )
 324 *
 325 *        Loop through eigenvalues of block nblk.
 326 *
 327    60    CONTINUE
 328          JBLK = 0
 329          DO 170 J = J1, M
 330             IF( IBLOCK( J ).NE.NBLK ) THEN
 331                J1 = J
 332                GO TO 180
 333             END IF
 334             JBLK = JBLK + 1
 335             XJ = W( J )
 336 *
 337 *           Skip all the work if the block size is one.
 338 *
 339             IF( BLKSIZ.EQ.1 ) THEN
 340                WORK( INDRV1+1 ) = ONE
 341                GO TO 140
 342             END IF
 343 *
 344 *           If eigenvalues j and j-1 are too close, add a relatively
 345 *           small perturbation.
 346 *
 347             IF( JBLK.GT.1 ) THEN
 348                EPS1 = ABS( EPS*XJ )
 349                PERTOL = TEN*EPS1
 350                SEP = XJ - XJM
 351                IF( SEP.LT.PERTOL )
 352      $            XJ = XJM + PERTOL
 353             END IF
 354 *
 355             ITS = 0
 356             NRMCHK = 0
 357 *
 358 *           Get random starting vector.
 359 *
 360             CALL DLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) )
 361 *
 362 *           Copy the matrix T so it won't be destroyed in factorization.
 363 *
 364             CALL DCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 )
 365             CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 )
 366             CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV3+1 ), 1 )
 367 *
 368 *           Compute LU factors with partial pivoting  ( PT = LU )
 369 *
 370             TOL = ZERO
 371             CALL DLAGTF( BLKSIZ, WORK( INDRV4+1 ), XJ, WORK( INDRV2+2 ),
 372      $                   WORK( INDRV3+1 ), TOL, WORK( INDRV5+1 ), IWORK,
 373      $                   IINFO )
 374 *
 375 *           Update iteration count.
 376 *
 377    70       CONTINUE
 378             ITS = ITS + 1
 379             IF( ITS.GT.MAXITS )
 380      $         GO TO 120
 381 *
 382 *           Normalize and scale the righthand side vector Pb.
 383 *
 384             JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
 385             SCL = BLKSIZ*ONENRM*MAX( EPS,
 386      $            ABS( WORK( INDRV4+BLKSIZ ) ) ) /
 387      $            ABS( WORK( INDRV1+JMAX ) )
 388             CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
 389 *
 390 *           Solve the system LU = Pb.
 391 *
 392             CALL DLAGTS( -1, BLKSIZ, WORK( INDRV4+1 ), WORK( INDRV2+2 ),
 393      $                   WORK( INDRV3+1 ), WORK( INDRV5+1 ), IWORK,
 394      $                   WORK( INDRV1+1 ), TOL, IINFO )
 395 *
 396 *           Reorthogonalize by modified Gram-Schmidt if eigenvalues are
 397 *           close enough.
 398 *
 399             IF( JBLK.EQ.1 )
 400      $         GO TO 110
 401             IF( ABS( XJ-XJM ).GT.ORTOL )
 402      $         GPIND = J
 403             IF( GPIND.NE.J ) THEN
 404                DO 100 I = GPIND, J - 1
 405                   ZTR = ZERO
 406                   DO 80 JR = 1, BLKSIZ
 407                      ZTR = ZTR + WORK( INDRV1+JR )*
 408      $                     DBLE( Z( B1-1+JR, I ) )
 409    80             CONTINUE
 410                   DO 90 JR = 1, BLKSIZ
 411                      WORK( INDRV1+JR ) = WORK( INDRV1+JR ) -
 412      $                                   ZTR*DBLE( Z( B1-1+JR, I ) )
 413    90             CONTINUE
 414   100          CONTINUE
 415             END IF
 416 *
 417 *           Check the infinity norm of the iterate.
 418 *
 419   110       CONTINUE
 420             JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
 421             NRM = ABS( WORK( INDRV1+JMAX ) )
 422 *
 423 *           Continue for additional iterations after norm reaches
 424 *           stopping criterion.
 425 *
 426             IF( NRM.LT.DTPCRT )
 427      $         GO TO 70
 428             NRMCHK = NRMCHK + 1
 429             IF( NRMCHK.LT.EXTRA+1 )
 430      $         GO TO 70
 431 *
 432             GO TO 130
 433 *
 434 *           If stopping criterion was not satisfied, update info and
 435 *           store eigenvector number in array ifail.
 436 *
 437   120       CONTINUE
 438             INFO = INFO + 1
 439             IFAIL( INFO ) = J
 440 *
 441 *           Accept iterate as jth eigenvector.
 442 *
 443   130       CONTINUE
 444             SCL = ONE / DNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 )
 445             JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
 446             IF( WORK( INDRV1+JMAX ).LT.ZERO )
 447      $         SCL = -SCL
 448             CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
 449   140       CONTINUE
 450             DO 150 I = 1, N
 451                Z( I, J ) = CZERO
 452   150       CONTINUE
 453             DO 160 I = 1, BLKSIZ
 454                Z( B1+I-1, J ) = DCMPLX( WORK( INDRV1+I ), ZERO )
 455   160       CONTINUE
 456 *
 457 *           Save the shift to check eigenvalue spacing at next
 458 *           iteration.
 459 *
 460             XJM = XJ
 461 *
 462   170    CONTINUE
 463   180 CONTINUE
 464 *
 465       RETURN
 466 *
 467 *     End of ZSTEIN
 468 *
 469       END