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