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