1      SUBROUTINE DSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
2     $                   IWORK, IFAIL, INFO )
3*
4*  -- LAPACK routine (instrumented to count operations, version 3.0) --
5*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6*     Courant Institute, Argonne National Lab, and Rice University
7*     September 30, 1994
8*
9*     .. Scalar Arguments ..
10      INTEGER            INFO, LDZ, M, N
11*     ..
12*     .. Array Arguments ..
13      INTEGER            IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
14     $                   IWORK( * )
15      DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
16*     ..
17*     Common block to return operation count and iteration count
18*     ITCNT is initialized to 0, OPS is only incremented
19*     .. Common blocks ..
20      COMMON             / LATIME / OPS, ITCNT
21*     ..
22*     .. Scalars in Common ..
23      DOUBLE PRECISION   ITCNT, OPS
24*     ..
25*
26*  Purpose
27*  =======
28*
29*  DSTEIN computes the eigenvectors of a real symmetric tridiagonal
30*  matrix T corresponding to specified eigenvalues, using inverse
31*  iteration.
32*
33*  The maximum number of iterations allowed for each eigenvector is
34*  specified by an internal parameter MAXITS (currently set to 5).
35*
36*  Arguments
37*  =========
38*
39*  N       (input) INTEGER
40*          The order of the matrix.  N >= 0.
41*
42*  D       (input) DOUBLE PRECISION array, dimension (N)
43*          The n diagonal elements of the tridiagonal matrix T.
44*
45*  E       (input) DOUBLE PRECISION array, dimension (N)
46*          The (n-1) subdiagonal elements of the tridiagonal matrix
47*          T, in elements 1 to N-1.  E(N) need not be set.
48*
49*  M       (input) INTEGER
50*          The number of eigenvectors to be found.  0 <= M <= N.
51*
52*  W       (input) DOUBLE PRECISION array, dimension (N)
53*          The first M elements of W contain the eigenvalues for
54*          which eigenvectors are to be computed.  The eigenvalues
55*          should be grouped by split-off block and ordered from
56*          smallest to largest within the block.  ( The output array
57*          W from DSTEBZ with ORDER = 'B' is expected here. )
58*
59*  IBLOCK  (input) INTEGER array, dimension (N)
60*          The submatrix indices associated with the corresponding
61*          eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
62*          the first submatrix from the top, =2 if W(i) belongs to
63*          the second submatrix, etc.  ( The output array IBLOCK
64*          from DSTEBZ is expected here. )
65*
66*  ISPLIT  (input) INTEGER array, dimension (N)
67*          The splitting points, at which T breaks up into submatrices.
68*          The first submatrix consists of rows/columns 1 to
69*          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
70*          through ISPLIT( 2 ), etc.
71*          ( The output array ISPLIT from DSTEBZ is expected here. )
72*
73*  Z       (output) DOUBLE PRECISION array, dimension (LDZ, M)
74*          The computed eigenvectors.  The eigenvector associated
75*          with the eigenvalue W(i) is stored in the i-th column of
76*          Z.  Any vector which fails to converge is set to its current
77*          iterate after MAXITS iterations.
78*
79*  LDZ     (input) INTEGER
80*          The leading dimension of the array Z.  LDZ >= max(1,N).
81*
82*  WORK    (workspace) DOUBLE PRECISION array, dimension (5*N)
83*
84*  IWORK   (workspace) INTEGER array, dimension (N)
85*
86*  IFAIL   (output) INTEGER array, dimension (M)
87*          On normal exit, all elements of IFAIL are zero.
88*          If one or more eigenvectors fail to converge after
89*          MAXITS iterations, then their indices are stored in
90*          array IFAIL.
91*
92*  INFO    (output) INTEGER
93*          = 0: successful exit.
94*          < 0: if INFO = -i, the i-th argument had an illegal value
95*          > 0: if INFO = i, then i eigenvectors failed to converge
96*               in MAXITS iterations.  Their indices are stored in
97*               array IFAIL.
98*
99*  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      DOUBLE PRECISION   ZERO, ONE, TEN, ODM3, ODM1
113      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1,
114     $                   ODM3 = 1.0D-3, ODM1 = 1.0D-1 )
115      INTEGER            MAXITS, EXTRA
116      PARAMETER          ( MAXITS = 5, EXTRA = 2 )
117*     ..
118*     .. Local Scalars ..
119      INTEGER            B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
120     $                   INDRV2, INDRV3, INDRV4, INDRV5, ITS, J, J1,
121     $                   JBLK, JMAX, NBLK, NRMCHK
122      DOUBLE PRECISION   DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
123     $                   SCL, SEP, TOL, XJ, XJM, ZTR
124*     ..
125*     .. Local Arrays ..
126      INTEGER            ISEED( 4 )
127*     ..
128*     .. External Functions ..
129      INTEGER            IDAMAX
130      DOUBLE PRECISION   DASUM, DDOT, DLAMCH, DNRM2
131      EXTERNAL           IDAMAX, DASUM, DDOT, DLAMCH, DNRM2
132*     ..
133*     .. External Subroutines ..
134      EXTERNAL           DAXPY, DCOPY, DLAGTF, DLAGTS, DLARNV, DSCAL,
135     $                   XERBLA
136*     ..
137*     .. Intrinsic Functions ..
138      INTRINSIC          ABS, MAX, SQRT
139*     ..
140*     .. Executable Statements ..
141*
142*     Test the input parameters.
143*
144      INFO = 0
145      DO 10 I = 1, M
146         IFAIL( I ) = 0
147   10 CONTINUE
148*
149      IF( N.LT.0 ) THEN
150         INFO = -1
151      ELSE IF( M.LT.0 .OR. M.GT.N ) THEN
152         INFO = -4
153      ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN
154         INFO = -9
155      ELSE
156         DO 20 J = 2, M
157            IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN
158               INFO = -6
159               GO TO 30
160            END IF
161            IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) )
162     $           THEN
163               INFO = -5
164               GO TO 30
165            END IF
166   20    CONTINUE
167   30    CONTINUE
168      END IF
169*
170      IF( INFO.NE.0 ) THEN
171         CALL XERBLA( 'DSTEIN', -INFO )
172         RETURN
173      END IF
174*
175*     Initialize iteration count.
176*
177      ITCNT = 0
178*
179*     Quick return if possible
180*
181      IF( N.EQ.0 .OR. M.EQ.0 ) THEN
182         RETURN
183      ELSE IF( N.EQ.1 ) THEN
184         Z( 1, 1 ) = ONE
185         RETURN
186      END IF
187*
188*     Get machine constants.
189*
190      EPS = DLAMCH( 'Precision' )
191*
192*     Initialize seed for random number generator DLARNV.
193*
194      DO 40 I = 1, 4
195         ISEED( I ) = 1
196   40 CONTINUE
197*
198*     Initialize pointers.
199*
200      INDRV1 = 0
201      INDRV2 = INDRV1 + N
202      INDRV3 = INDRV2 + N
203      INDRV4 = INDRV3 + N
204      INDRV5 = INDRV4 + N
205*
206*     Compute eigenvectors of matrix blocks.
207*
208      J1 = 1
209      DO 160 NBLK = 1, IBLOCK( M )
210*
211*        Find starting and ending indices of block nblk.
212*
213         IF( NBLK.EQ.1 ) THEN
214            B1 = 1
215         ELSE
216            B1 = ISPLIT( NBLK-1 ) + 1
217         END IF
218         BN = ISPLIT( NBLK )
219         BLKSIZ = BN - B1 + 1
220         IF( BLKSIZ.EQ.1 )
221     $      GO TO 60
222         GPIND = B1
223*
224*        Compute reorthogonalization criterion and stopping criterion.
225*
226         ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) )
227         ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) )
228         DO 50 I = B1 + 1, BN - 1
229            ONENRM = MAX( ONENRM, ABS( D( I ) )+ABS( E( I-1 ) )+
230     $               ABS( E( I ) ) )
231   50    CONTINUE
232         ORTOL = ODM3*ONENRM
233*
234         DTPCRT = SQRT( ODM1 / BLKSIZ )
235*
236*        Increment opcount for computing criteria.
237*
238         OPS = OPS + ( BN-B1 )*2 + 3
239*
240*        Loop through eigenvalues of block nblk.
241*
242   60    CONTINUE
243         JBLK = 0
244         DO 150 J = J1, M
245            IF( IBLOCK( J ).NE.NBLK ) THEN
246               J1 = J
247               GO TO 160
248            END IF
249            JBLK = JBLK + 1
250            XJ = W( J )
251*
252*           Skip all the work if the block size is one.
253*
254            IF( BLKSIZ.EQ.1 ) THEN
255               WORK( INDRV1+1 ) = ONE
256               GO TO 120
257            END IF
258*
259*           If eigenvalues j and j-1 are too close, add a relatively
260*           small perturbation.
261*
262            IF( JBLK.GT.1 ) THEN
263               EPS1 = ABS( EPS*XJ )
264               PERTOL = TEN*EPS1
265               SEP = XJ - XJM
266               IF( SEP.LT.PERTOL )
267     $            XJ = XJM + PERTOL
268            END IF
269*
270            ITS = 0
271            NRMCHK = 0
272*
273*           Get random starting vector.
274*
275            CALL DLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) )
276*
277*           Increment opcount for getting random starting vector.
278*           ( DLARND(2,.) requires 9 flops. )
279*
280            OPS = OPS + BLKSIZ*9
281*
282*           Copy the matrix T so it won't be destroyed in factorization.
283*
284            CALL DCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 )
285            CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 )
286            CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV3+1 ), 1 )
287*
288*           Compute LU factors with partial pivoting  ( PT = LU )
289*
290            TOL = ZERO
291            CALL DLAGTF( BLKSIZ, WORK( INDRV4+1 ), XJ, WORK( INDRV2+2 ),
292     $                   WORK( INDRV3+1 ), TOL, WORK( INDRV5+1 ), IWORK,
293     $                   IINFO )
294*
295*           Increment opcount for computing LU factors.
296*           ( DLAGTF(BLKSIZ,...) requires about 8*BLKSIZ flops. )
297*
298            OPS = OPS + 8*BLKSIZ
299*
300*           Update iteration count.
301*
302   70       CONTINUE
303            ITS = ITS + 1
304            IF( ITS.GT.MAXITS )
305     $         GO TO 100
306*
307*           Normalize and scale the righthand side vector Pb.
308*
309            SCL = BLKSIZ*ONENRM*MAX( EPS,
310     $            ABS( WORK( INDRV4+BLKSIZ ) ) ) /
311     $            DASUM( BLKSIZ, WORK( INDRV1+1 ), 1 )
312            CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
313*
314*           Solve the system LU = Pb.
315*
316            CALL DLAGTS( -1, BLKSIZ, WORK( INDRV4+1 ), WORK( INDRV2+2 ),
317     $                   WORK( INDRV3+1 ), WORK( INDRV5+1 ), IWORK,
318     $                   WORK( INDRV1+1 ), TOL, IINFO )
319*
320*           Increment opcount for scaling and solving linear system.
321*           ( DLAGTS(-1,BLKSIZ,...) requires about 8*BLKSIZ flops. )
322*
323            OPS = OPS + 3 + 10*BLKSIZ
324*
325*           Reorthogonalize by modified Gram-Schmidt if eigenvalues are
326*           close enough.
327*
328            IF( JBLK.EQ.1 )
329     $         GO TO 90
330            IF( ABS( XJ-XJM ).GT.ORTOL )
331     $         GPIND = J
332            IF( GPIND.NE.J ) THEN
333               DO 80 I = GPIND, J - 1
334                  ZTR = -DDOT( BLKSIZ, WORK( INDRV1+1 ), 1, Z( B1, I ),
335     $                  1 )
336                  CALL DAXPY( BLKSIZ, ZTR, Z( B1, I ), 1,
337     $                        WORK( INDRV1+1 ), 1 )
338   80          CONTINUE
339*
340*              Increment opcount for reorthogonalizing.
341*
342               OPS = OPS + ( J-GPIND )*BLKSIZ*4
343*
344            END IF
345*
346*           Check the infinity norm of the iterate.
347*
348   90       CONTINUE
349            JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
350            NRM = ABS( WORK( INDRV1+JMAX ) )
351*
352*           Continue for additional iterations after norm reaches
353*           stopping criterion.
354*
355            IF( NRM.LT.DTPCRT )
356     $         GO TO 70
357            NRMCHK = NRMCHK + 1
358            IF( NRMCHK.LT.EXTRA+1 )
359     $         GO TO 70
360*
361            GO TO 110
362*
363*           If stopping criterion was not satisfied, update info and
364*           store eigenvector number in array ifail.
365*
366  100       CONTINUE
367            INFO = INFO + 1
368            IFAIL( INFO ) = J
369*
370*           Accept iterate as jth eigenvector.
371*
372  110       CONTINUE
373            SCL = ONE / DNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 )
374            JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
375            IF( WORK( INDRV1+JMAX ).LT.ZERO )
376     $         SCL = -SCL
377            CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
378*
379*           Increment opcount for scaling.
380*
381            OPS = OPS + 3*BLKSIZ
382*
383  120       CONTINUE
384            DO 130 I = 1, N
385               Z( I, J ) = ZERO
386  130       CONTINUE
387            DO 140 I = 1, BLKSIZ
388               Z( B1+I-1, J ) = WORK( INDRV1+I )
389  140       CONTINUE
390*
391*           Save the shift to check eigenvalue spacing at next
392*           iteration.
393*
394            XJM = XJ
395*
396  150    CONTINUE
397  160 CONTINUE
398*
399      RETURN
400*
401*     End of DSTEIN
402*
403      END
404