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*> \ingroup doubleOTHERcomputational 170* 171* ===================================================================== 172 SUBROUTINE DSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, 173 $ IWORK, IFAIL, INFO ) 174* 175* -- LAPACK computational routine -- 176* -- LAPACK is a software package provided by Univ. of Tennessee, -- 177* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 178* 179* .. Scalar Arguments .. 180 INTEGER INFO, LDZ, M, N 181* .. 182* .. Array Arguments .. 183 INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ), 184 $ IWORK( * ) 185 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * ) 186* .. 187* 188* ===================================================================== 189* 190* .. Parameters .. 191 DOUBLE PRECISION ZERO, ONE, TEN, ODM3, ODM1 192 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1, 193 $ ODM3 = 1.0D-3, ODM1 = 1.0D-1 ) 194 INTEGER MAXITS, EXTRA 195 PARAMETER ( MAXITS = 5, EXTRA = 2 ) 196* .. 197* .. Local Scalars .. 198 INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1, 199 $ INDRV2, INDRV3, INDRV4, INDRV5, ITS, J, J1, 200 $ JBLK, JMAX, NBLK, NRMCHK 201 DOUBLE PRECISION DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL, 202 $ SCL, SEP, TOL, XJ, XJM, ZTR 203* .. 204* .. Local Arrays .. 205 INTEGER ISEED( 4 ) 206* .. 207* .. External Functions .. 208 INTEGER IDAMAX 209 DOUBLE PRECISION DDOT, DLAMCH, DNRM2 210 EXTERNAL IDAMAX, DDOT, DLAMCH, DNRM2 211* .. 212* .. External Subroutines .. 213 EXTERNAL DAXPY, DCOPY, DLAGTF, DLAGTS, DLARNV, DSCAL, 214 $ XERBLA 215* .. 216* .. Intrinsic Functions .. 217 INTRINSIC ABS, MAX, SQRT 218* .. 219* .. Executable Statements .. 220* 221* Test the input parameters. 222* 223 INFO = 0 224 DO 10 I = 1, M 225 IFAIL( I ) = 0 226 10 CONTINUE 227* 228 IF( N.LT.0 ) THEN 229 INFO = -1 230 ELSE IF( M.LT.0 .OR. M.GT.N ) THEN 231 INFO = -4 232 ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN 233 INFO = -9 234 ELSE 235 DO 20 J = 2, M 236 IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN 237 INFO = -6 238 GO TO 30 239 END IF 240 IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) ) 241 $ THEN 242 INFO = -5 243 GO TO 30 244 END IF 245 20 CONTINUE 246 30 CONTINUE 247 END IF 248* 249 IF( INFO.NE.0 ) THEN 250 CALL XERBLA( 'DSTEIN', -INFO ) 251 RETURN 252 END IF 253* 254* Quick return if possible 255* 256 IF( N.EQ.0 .OR. M.EQ.0 ) THEN 257 RETURN 258 ELSE IF( N.EQ.1 ) THEN 259 Z( 1, 1 ) = ONE 260 RETURN 261 END IF 262* 263* Get machine constants. 264* 265 EPS = DLAMCH( 'Precision' ) 266* 267* Initialize seed for random number generator DLARNV. 268* 269 DO 40 I = 1, 4 270 ISEED( I ) = 1 271 40 CONTINUE 272* 273* Initialize pointers. 274* 275 INDRV1 = 0 276 INDRV2 = INDRV1 + N 277 INDRV3 = INDRV2 + N 278 INDRV4 = INDRV3 + N 279 INDRV5 = INDRV4 + N 280* 281* Compute eigenvectors of matrix blocks. 282* 283 J1 = 1 284 DO 160 NBLK = 1, IBLOCK( M ) 285* 286* Find starting and ending indices of block nblk. 287* 288 IF( NBLK.EQ.1 ) THEN 289 B1 = 1 290 ELSE 291 B1 = ISPLIT( NBLK-1 ) + 1 292 END IF 293 BN = ISPLIT( NBLK ) 294 BLKSIZ = BN - B1 + 1 295 IF( BLKSIZ.EQ.1 ) 296 $ GO TO 60 297 GPIND = J1 298* 299* Compute reorthogonalization criterion and stopping criterion. 300* 301 ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) ) 302 ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) ) 303 DO 50 I = B1 + 1, BN - 1 304 ONENRM = MAX( ONENRM, ABS( D( I ) )+ABS( E( I-1 ) )+ 305 $ ABS( E( I ) ) ) 306 50 CONTINUE 307 ORTOL = ODM3*ONENRM 308* 309 DTPCRT = SQRT( ODM1 / BLKSIZ ) 310* 311* Loop through eigenvalues of block nblk. 312* 313 60 CONTINUE 314 JBLK = 0 315 DO 150 J = J1, M 316 IF( IBLOCK( J ).NE.NBLK ) THEN 317 J1 = J 318 GO TO 160 319 END IF 320 JBLK = JBLK + 1 321 XJ = W( J ) 322* 323* Skip all the work if the block size is one. 324* 325 IF( BLKSIZ.EQ.1 ) THEN 326 WORK( INDRV1+1 ) = ONE 327 GO TO 120 328 END IF 329* 330* If eigenvalues j and j-1 are too close, add a relatively 331* small perturbation. 332* 333 IF( JBLK.GT.1 ) THEN 334 EPS1 = ABS( EPS*XJ ) 335 PERTOL = TEN*EPS1 336 SEP = XJ - XJM 337 IF( SEP.LT.PERTOL ) 338 $ XJ = XJM + PERTOL 339 END IF 340* 341 ITS = 0 342 NRMCHK = 0 343* 344* Get random starting vector. 345* 346 CALL DLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) ) 347* 348* Copy the matrix T so it won't be destroyed in factorization. 349* 350 CALL DCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 ) 351 CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 ) 352 CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV3+1 ), 1 ) 353* 354* Compute LU factors with partial pivoting ( PT = LU ) 355* 356 TOL = ZERO 357 CALL DLAGTF( BLKSIZ, WORK( INDRV4+1 ), XJ, WORK( INDRV2+2 ), 358 $ WORK( INDRV3+1 ), TOL, WORK( INDRV5+1 ), IWORK, 359 $ IINFO ) 360* 361* Update iteration count. 362* 363 70 CONTINUE 364 ITS = ITS + 1 365 IF( ITS.GT.MAXITS ) 366 $ GO TO 100 367* 368* Normalize and scale the righthand side vector Pb. 369* 370 JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 ) 371 SCL = BLKSIZ*ONENRM*MAX( EPS, 372 $ ABS( WORK( INDRV4+BLKSIZ ) ) ) / 373 $ ABS( WORK( INDRV1+JMAX ) ) 374 CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 ) 375* 376* Solve the system LU = Pb. 377* 378 CALL DLAGTS( -1, BLKSIZ, WORK( INDRV4+1 ), WORK( INDRV2+2 ), 379 $ WORK( INDRV3+1 ), WORK( INDRV5+1 ), IWORK, 380 $ WORK( INDRV1+1 ), TOL, IINFO ) 381* 382* Reorthogonalize by modified Gram-Schmidt if eigenvalues are 383* close enough. 384* 385 IF( JBLK.EQ.1 ) 386 $ GO TO 90 387 IF( ABS( XJ-XJM ).GT.ORTOL ) 388 $ GPIND = J 389 IF( GPIND.NE.J ) THEN 390 DO 80 I = GPIND, J - 1 391 ZTR = -DDOT( BLKSIZ, WORK( INDRV1+1 ), 1, Z( B1, I ), 392 $ 1 ) 393 CALL DAXPY( BLKSIZ, ZTR, Z( B1, I ), 1, 394 $ WORK( INDRV1+1 ), 1 ) 395 80 CONTINUE 396 END IF 397* 398* Check the infinity norm of the iterate. 399* 400 90 CONTINUE 401 JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 ) 402 NRM = ABS( WORK( INDRV1+JMAX ) ) 403* 404* Continue for additional iterations after norm reaches 405* stopping criterion. 406* 407 IF( NRM.LT.DTPCRT ) 408 $ GO TO 70 409 NRMCHK = NRMCHK + 1 410 IF( NRMCHK.LT.EXTRA+1 ) 411 $ GO TO 70 412* 413 GO TO 110 414* 415* If stopping criterion was not satisfied, update info and 416* store eigenvector number in array ifail. 417* 418 100 CONTINUE 419 INFO = INFO + 1 420 IFAIL( INFO ) = J 421* 422* Accept iterate as jth eigenvector. 423* 424 110 CONTINUE 425 SCL = ONE / DNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 ) 426 JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 ) 427 IF( WORK( INDRV1+JMAX ).LT.ZERO ) 428 $ SCL = -SCL 429 CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 ) 430 120 CONTINUE 431 DO 130 I = 1, N 432 Z( I, J ) = ZERO 433 130 CONTINUE 434 DO 140 I = 1, BLKSIZ 435 Z( B1+I-1, J ) = WORK( INDRV1+I ) 436 140 CONTINUE 437* 438* Save the shift to check eigenvalue spacing at next 439* iteration. 440* 441 XJM = XJ 442* 443 150 CONTINUE 444 160 CONTINUE 445* 446 RETURN 447* 448* End of DSTEIN 449* 450 END 451