001:       SUBROUTINE DSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
002:      $                   IWORK, IFAIL, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            INFO, LDZ, M, N
011: *     ..
012: *     .. Array Arguments ..
013:       INTEGER            IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
014:      $                   IWORK( * )
015:       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  DSTEIN computes the eigenvectors of a real symmetric tridiagonal
022: *  matrix T corresponding to specified eigenvalues, using inverse
023: *  iteration.
024: *
025: *  The maximum number of iterations allowed for each eigenvector is
026: *  specified by an internal parameter MAXITS (currently set to 5).
027: *
028: *  Arguments
029: *  =========
030: *
031: *  N       (input) INTEGER
032: *          The order of the matrix.  N >= 0.
033: *
034: *  D       (input) DOUBLE PRECISION array, dimension (N)
035: *          The n diagonal elements of the tridiagonal matrix T.
036: *
037: *  E       (input) DOUBLE PRECISION array, dimension (N-1)
038: *          The (n-1) subdiagonal elements of the tridiagonal matrix
039: *          T, in elements 1 to N-1.
040: *
041: *  M       (input) INTEGER
042: *          The number of eigenvectors to be found.  0 <= M <= N.
043: *
044: *  W       (input) DOUBLE PRECISION array, dimension (N)
045: *          The first M elements of W contain the eigenvalues for
046: *          which eigenvectors are to be computed.  The eigenvalues
047: *          should be grouped by split-off block and ordered from
048: *          smallest to largest within the block.  ( The output array
049: *          W from DSTEBZ with ORDER = 'B' is expected here. )
050: *
051: *  IBLOCK  (input) INTEGER array, dimension (N)
052: *          The submatrix indices associated with the corresponding
053: *          eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
054: *          the first submatrix from the top, =2 if W(i) belongs to
055: *          the second submatrix, etc.  ( The output array IBLOCK
056: *          from DSTEBZ is expected here. )
057: *
058: *  ISPLIT  (input) INTEGER array, dimension (N)
059: *          The splitting points, at which T breaks up into submatrices.
060: *          The first submatrix consists of rows/columns 1 to
061: *          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
062: *          through ISPLIT( 2 ), etc.
063: *          ( The output array ISPLIT from DSTEBZ is expected here. )
064: *
065: *  Z       (output) DOUBLE PRECISION array, dimension (LDZ, M)
066: *          The computed eigenvectors.  The eigenvector associated
067: *          with the eigenvalue W(i) is stored in the i-th column of
068: *          Z.  Any vector which fails to converge is set to its current
069: *          iterate after MAXITS iterations.
070: *
071: *  LDZ     (input) INTEGER
072: *          The leading dimension of the array Z.  LDZ >= max(1,N).
073: *
074: *  WORK    (workspace) DOUBLE PRECISION array, dimension (5*N)
075: *
076: *  IWORK   (workspace) INTEGER array, dimension (N)
077: *
078: *  IFAIL   (output) INTEGER array, dimension (M)
079: *          On normal exit, all elements of IFAIL are zero.
080: *          If one or more eigenvectors fail to converge after
081: *          MAXITS iterations, then their indices are stored in
082: *          array IFAIL.
083: *
084: *  INFO    (output) INTEGER
085: *          = 0: successful exit.
086: *          < 0: if INFO = -i, the i-th argument had an illegal value
087: *          > 0: if INFO = i, then i eigenvectors failed to converge
088: *               in MAXITS iterations.  Their indices are stored in
089: *               array IFAIL.
090: *
091: *  Internal Parameters
092: *  ===================
093: *
094: *  MAXITS  INTEGER, default = 5
095: *          The maximum number of iterations performed.
096: *
097: *  EXTRA   INTEGER, default = 2
098: *          The number of iterations performed after norm growth
099: *          criterion is satisfied, should be at least 1.
100: *
101: *  =====================================================================
102: *
103: *     .. Parameters ..
104:       DOUBLE PRECISION   ZERO, ONE, TEN, ODM3, ODM1
105:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1,
106:      $                   ODM3 = 1.0D-3, ODM1 = 1.0D-1 )
107:       INTEGER            MAXITS, EXTRA
108:       PARAMETER          ( MAXITS = 5, EXTRA = 2 )
109: *     ..
110: *     .. Local Scalars ..
111:       INTEGER            B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
112:      $                   INDRV2, INDRV3, INDRV4, INDRV5, ITS, J, J1,
113:      $                   JBLK, JMAX, NBLK, NRMCHK
114:       DOUBLE PRECISION   DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
115:      $                   SCL, SEP, TOL, XJ, XJM, ZTR
116: *     ..
117: *     .. Local Arrays ..
118:       INTEGER            ISEED( 4 )
119: *     ..
120: *     .. External Functions ..
121:       INTEGER            IDAMAX
122:       DOUBLE PRECISION   DASUM, DDOT, DLAMCH, DNRM2
123:       EXTERNAL           IDAMAX, DASUM, DDOT, DLAMCH, DNRM2
124: *     ..
125: *     .. External Subroutines ..
126:       EXTERNAL           DAXPY, DCOPY, DLAGTF, DLAGTS, DLARNV, DSCAL,
127:      $                   XERBLA
128: *     ..
129: *     .. Intrinsic Functions ..
130:       INTRINSIC          ABS, MAX, SQRT
131: *     ..
132: *     .. Executable Statements ..
133: *
134: *     Test the input parameters.
135: *
136:       INFO = 0
137:       DO 10 I = 1, M
138:          IFAIL( I ) = 0
139:    10 CONTINUE
140: *
141:       IF( N.LT.0 ) THEN
142:          INFO = -1
143:       ELSE IF( M.LT.0 .OR. M.GT.N ) THEN
144:          INFO = -4
145:       ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN
146:          INFO = -9
147:       ELSE
148:          DO 20 J = 2, M
149:             IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN
150:                INFO = -6
151:                GO TO 30
152:             END IF
153:             IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) )
154:      $           THEN
155:                INFO = -5
156:                GO TO 30
157:             END IF
158:    20    CONTINUE
159:    30    CONTINUE
160:       END IF
161: *
162:       IF( INFO.NE.0 ) THEN
163:          CALL XERBLA( 'DSTEIN', -INFO )
164:          RETURN
165:       END IF
166: *
167: *     Quick return if possible
168: *
169:       IF( N.EQ.0 .OR. M.EQ.0 ) THEN
170:          RETURN
171:       ELSE IF( N.EQ.1 ) THEN
172:          Z( 1, 1 ) = ONE
173:          RETURN
174:       END IF
175: *
176: *     Get machine constants.
177: *
178:       EPS = DLAMCH( 'Precision' )
179: *
180: *     Initialize seed for random number generator DLARNV.
181: *
182:       DO 40 I = 1, 4
183:          ISEED( I ) = 1
184:    40 CONTINUE
185: *
186: *     Initialize pointers.
187: *
188:       INDRV1 = 0
189:       INDRV2 = INDRV1 + N
190:       INDRV3 = INDRV2 + N
191:       INDRV4 = INDRV3 + N
192:       INDRV5 = INDRV4 + N
193: *
194: *     Compute eigenvectors of matrix blocks.
195: *
196:       J1 = 1
197:       DO 160 NBLK = 1, IBLOCK( M )
198: *
199: *        Find starting and ending indices of block nblk.
200: *
201:          IF( NBLK.EQ.1 ) THEN
202:             B1 = 1
203:          ELSE
204:             B1 = ISPLIT( NBLK-1 ) + 1
205:          END IF
206:          BN = ISPLIT( NBLK )
207:          BLKSIZ = BN - B1 + 1
208:          IF( BLKSIZ.EQ.1 )
209:      $      GO TO 60
210:          GPIND = B1
211: *
212: *        Compute reorthogonalization criterion and stopping criterion.
213: *
214:          ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) )
215:          ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) )
216:          DO 50 I = B1 + 1, BN - 1
217:             ONENRM = MAX( ONENRM, ABS( D( I ) )+ABS( E( I-1 ) )+
218:      $               ABS( E( I ) ) )
219:    50    CONTINUE
220:          ORTOL = ODM3*ONENRM
221: *
222:          DTPCRT = SQRT( ODM1 / BLKSIZ )
223: *
224: *        Loop through eigenvalues of block nblk.
225: *
226:    60    CONTINUE
227:          JBLK = 0
228:          DO 150 J = J1, M
229:             IF( IBLOCK( J ).NE.NBLK ) THEN
230:                J1 = J
231:                GO TO 160
232:             END IF
233:             JBLK = JBLK + 1
234:             XJ = W( J )
235: *
236: *           Skip all the work if the block size is one.
237: *
238:             IF( BLKSIZ.EQ.1 ) THEN
239:                WORK( INDRV1+1 ) = ONE
240:                GO TO 120
241:             END IF
242: *
243: *           If eigenvalues j and j-1 are too close, add a relatively
244: *           small perturbation.
245: *
246:             IF( JBLK.GT.1 ) THEN
247:                EPS1 = ABS( EPS*XJ )
248:                PERTOL = TEN*EPS1
249:                SEP = XJ - XJM
250:                IF( SEP.LT.PERTOL )
251:      $            XJ = XJM + PERTOL
252:             END IF
253: *
254:             ITS = 0
255:             NRMCHK = 0
256: *
257: *           Get random starting vector.
258: *
259:             CALL DLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) )
260: *
261: *           Copy the matrix T so it won't be destroyed in factorization.
262: *
263:             CALL DCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 )
264:             CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 )
265:             CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV3+1 ), 1 )
266: *
267: *           Compute LU factors with partial pivoting  ( PT = LU )
268: *
269:             TOL = ZERO
270:             CALL DLAGTF( BLKSIZ, WORK( INDRV4+1 ), XJ, WORK( INDRV2+2 ),
271:      $                   WORK( INDRV3+1 ), TOL, WORK( INDRV5+1 ), IWORK,
272:      $                   IINFO )
273: *
274: *           Update iteration count.
275: *
276:    70       CONTINUE
277:             ITS = ITS + 1
278:             IF( ITS.GT.MAXITS )
279:      $         GO TO 100
280: *
281: *           Normalize and scale the righthand side vector Pb.
282: *
283:             SCL = BLKSIZ*ONENRM*MAX( EPS,
284:      $            ABS( WORK( INDRV4+BLKSIZ ) ) ) /
285:      $            DASUM( BLKSIZ, WORK( INDRV1+1 ), 1 )
286:             CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
287: *
288: *           Solve the system LU = Pb.
289: *
290:             CALL DLAGTS( -1, BLKSIZ, WORK( INDRV4+1 ), WORK( INDRV2+2 ),
291:      $                   WORK( INDRV3+1 ), WORK( INDRV5+1 ), IWORK,
292:      $                   WORK( INDRV1+1 ), TOL, IINFO )
293: *
294: *           Reorthogonalize by modified Gram-Schmidt if eigenvalues are
295: *           close enough.
296: *
297:             IF( JBLK.EQ.1 )
298:      $         GO TO 90
299:             IF( ABS( XJ-XJM ).GT.ORTOL )
300:      $         GPIND = J
301:             IF( GPIND.NE.J ) THEN
302:                DO 80 I = GPIND, J - 1
303:                   ZTR = -DDOT( BLKSIZ, WORK( INDRV1+1 ), 1, Z( B1, I ),
304:      $                  1 )
305:                   CALL DAXPY( BLKSIZ, ZTR, Z( B1, I ), 1,
306:      $                        WORK( INDRV1+1 ), 1 )
307:    80          CONTINUE
308:             END IF
309: *
310: *           Check the infinity norm of the iterate.
311: *
312:    90       CONTINUE
313:             JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
314:             NRM = ABS( WORK( INDRV1+JMAX ) )
315: *
316: *           Continue for additional iterations after norm reaches
317: *           stopping criterion.
318: *
319:             IF( NRM.LT.DTPCRT )
320:      $         GO TO 70
321:             NRMCHK = NRMCHK + 1
322:             IF( NRMCHK.LT.EXTRA+1 )
323:      $         GO TO 70
324: *
325:             GO TO 110
326: *
327: *           If stopping criterion was not satisfied, update info and
328: *           store eigenvector number in array ifail.
329: *
330:   100       CONTINUE
331:             INFO = INFO + 1
332:             IFAIL( INFO ) = J
333: *
334: *           Accept iterate as jth eigenvector.
335: *
336:   110       CONTINUE
337:             SCL = ONE / DNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 )
338:             JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
339:             IF( WORK( INDRV1+JMAX ).LT.ZERO )
340:      $         SCL = -SCL
341:             CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
342:   120       CONTINUE
343:             DO 130 I = 1, N
344:                Z( I, J ) = ZERO
345:   130       CONTINUE
346:             DO 140 I = 1, BLKSIZ
347:                Z( B1+I-1, J ) = WORK( INDRV1+I )
348:   140       CONTINUE
349: *
350: *           Save the shift to check eigenvalue spacing at next
351: *           iteration.
352: *
353:             XJM = XJ
354: *
355:   150    CONTINUE
356:   160 CONTINUE
357: *
358:       RETURN
359: *
360: *     End of DSTEIN
361: *
362:       END
363: