SRC/dstein.f

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