001:       SUBROUTINE SSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
002:      $                   IWORK, IFAIL, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            INFO, LDZ, M, N
010: *     ..
011: *     .. Array Arguments ..
012:       INTEGER            IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
013:      $                   IWORK( * )
014:       REAL               D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  SSTEIN computes the eigenvectors of a real symmetric tridiagonal
021: *  matrix T corresponding to specified eigenvalues, using inverse
022: *  iteration.
023: *
024: *  The maximum number of iterations allowed for each eigenvector is
025: *  specified by an internal parameter MAXITS (currently set to 5).
026: *
027: *  Arguments
028: *  =========
029: *
030: *  N       (input) INTEGER
031: *          The order of the matrix.  N >= 0.
032: *
033: *  D       (input) REAL array, dimension (N)
034: *          The n diagonal elements of the tridiagonal matrix T.
035: *
036: *  E       (input) REAL array, dimension (N-1)
037: *          The (n-1) subdiagonal elements of the tridiagonal matrix
038: *          T, in elements 1 to N-1.
039: *
040: *  M       (input) INTEGER
041: *          The number of eigenvectors to be found.  0 <= M <= N.
042: *
043: *  W       (input) REAL array, dimension (N)
044: *          The first M elements of W contain the eigenvalues for
045: *          which eigenvectors are to be computed.  The eigenvalues
046: *          should be grouped by split-off block and ordered from
047: *          smallest to largest within the block.  ( The output array
048: *          W from SSTEBZ with ORDER = 'B' is expected here. )
049: *
050: *  IBLOCK  (input) INTEGER array, dimension (N)
051: *          The submatrix indices associated with the corresponding
052: *          eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
053: *          the first submatrix from the top, =2 if W(i) belongs to
054: *          the second submatrix, etc.  ( The output array IBLOCK
055: *          from SSTEBZ is expected here. )
056: *
057: *  ISPLIT  (input) INTEGER array, dimension (N)
058: *          The splitting points, at which T breaks up into submatrices.
059: *          The first submatrix consists of rows/columns 1 to
060: *          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
061: *          through ISPLIT( 2 ), etc.
062: *          ( The output array ISPLIT from SSTEBZ is expected here. )
063: *
064: *  Z       (output) REAL array, dimension (LDZ, M)
065: *          The computed eigenvectors.  The eigenvector associated
066: *          with the eigenvalue W(i) is stored in the i-th column of
067: *          Z.  Any vector which fails to converge is set to its current
068: *          iterate after MAXITS iterations.
069: *
070: *  LDZ     (input) INTEGER
071: *          The leading dimension of the array Z.  LDZ >= max(1,N).
072: *
073: *  WORK    (workspace) REAL array, dimension (5*N)
074: *
075: *  IWORK   (workspace) INTEGER array, dimension (N)
076: *
077: *  IFAIL   (output) INTEGER array, dimension (M)
078: *          On normal exit, all elements of IFAIL are zero.
079: *          If one or more eigenvectors fail to converge after
080: *          MAXITS iterations, then their indices are stored in
081: *          array IFAIL.
082: *
083: *  INFO    (output) INTEGER
084: *          = 0: successful exit.
085: *          < 0: if INFO = -i, the i-th argument had an illegal value
086: *          > 0: if INFO = i, then i eigenvectors failed to converge
087: *               in MAXITS iterations.  Their indices are stored in
088: *               array IFAIL.
089: *
090: *  Internal Parameters
091: *  ===================
092: *
093: *  MAXITS  INTEGER, default = 5
094: *          The maximum number of iterations performed.
095: *
096: *  EXTRA   INTEGER, default = 2
097: *          The number of iterations performed after norm growth
098: *          criterion is satisfied, should be at least 1.
099: *
100: *  =====================================================================
101: *
102: *     .. Parameters ..
103:       REAL               ZERO, ONE, TEN, ODM3, ODM1
104:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+1,
105:      $                   ODM3 = 1.0E-3, ODM1 = 1.0E-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:       REAL               CTR, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
114:      $                   SCL, SEP, STPCRT, TOL, XJ, XJM
115: *     ..
116: *     .. Local Arrays ..
117:       INTEGER            ISEED( 4 )
118: *     ..
119: *     .. External Functions ..
120:       INTEGER            ISAMAX
121:       REAL               SASUM, SDOT, SLAMCH, SNRM2
122:       EXTERNAL           ISAMAX, SASUM, SDOT, SLAMCH, SNRM2
123: *     ..
124: *     .. External Subroutines ..
125:       EXTERNAL           SAXPY, SCOPY, SLAGTF, SLAGTS, SLARNV, SSCAL,
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( 'SSTEIN', -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 = SLAMCH( 'Precision' )
178: *
179: *     Initialize seed for random number generator SLARNV.
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:          STPCRT = 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 SLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) )
259: *
260: *           Copy the matrix T so it won't be destroyed in factorization.
261: *
262:             CALL SCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 )
263:             CALL SCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 )
264:             CALL SCOPY( 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 SLAGTF( 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:      $            SASUM( BLKSIZ, WORK( INDRV1+1 ), 1 )
285:             CALL SSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
286: *
287: *           Solve the system LU = Pb.
288: *
289:             CALL SLAGTS( -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:                   CTR = -SDOT( BLKSIZ, WORK( INDRV1+1 ), 1, Z( B1, I ),
303:      $                  1 )
304:                   CALL SAXPY( BLKSIZ, CTR, 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 = ISAMAX( 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.STPCRT )
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 / SNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 )
337:             JMAX = ISAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
338:             IF( WORK( INDRV1+JMAX ).LT.ZERO )
339:      $         SCL = -SCL
340:             CALL SSCAL( 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 SSTEIN
360: *
361:       END
362: