src/LAPACK/zstein.f

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