src/LAPACK/dstein.f

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