1      SUBROUTINE SSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
2     $                   IWORK, IFAIL, INFO )
3*
4*  -- LAPACK routine (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      REAL               D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
16*     ..
17*
18*  Purpose
19*  =======
20*
21*  SSTEIN computes the eigenvectors of a real symmetric tridiagonal
22*  matrix T corresponding to specified eigenvalues, using inverse
23*  iteration.
24*
25*  The maximum number of iterations allowed for each eigenvector is
26*  specified by an internal parameter MAXITS (currently set to 5).
27*
28*  Arguments
29*  =========
30*
31*  N       (input) INTEGER
32*          The order of the matrix.  N >= 0.
33*
34*  D       (input) REAL array, dimension (N)
35*          The n diagonal elements of the tridiagonal matrix T.
36*
37*  E       (input) REAL array, dimension (N)
38*          The (n-1) subdiagonal elements of the tridiagonal matrix
39*          T, in elements 1 to N-1.  E(N) need not be set.
40*
41*  M       (input) INTEGER
42*          The number of eigenvectors to be found.  0 <= M <= N.
43*
44*  W       (input) REAL array, dimension (N)
45*          The first M elements of W contain the eigenvalues for
46*          which eigenvectors are to be computed.  The eigenvalues
47*          should be grouped by split-off block and ordered from
48*          smallest to largest within the block.  ( The output array
49*          W from SSTEBZ with ORDER = 'B' is expected here. )
50*
51*  IBLOCK  (input) INTEGER array, dimension (N)
52*          The submatrix indices associated with the corresponding
53*          eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
54*          the first submatrix from the top, =2 if W(i) belongs to
55*          the second submatrix, etc.  ( The output array IBLOCK
56*          from SSTEBZ is expected here. )
57*
58*  ISPLIT  (input) INTEGER array, dimension (N)
59*          The splitting points, at which T breaks up into submatrices.
60*          The first submatrix consists of rows/columns 1 to
61*          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
62*          through ISPLIT( 2 ), etc.
63*          ( The output array ISPLIT from SSTEBZ is expected here. )
64*
65*  Z       (output) REAL array, dimension (LDZ, M)
66*          The computed eigenvectors.  The eigenvector associated
67*          with the eigenvalue W(i) is stored in the i-th column of
68*          Z.  Any vector which fails to converge is set to its current
69*          iterate after MAXITS iterations.
70*
71*  LDZ     (input) INTEGER
72*          The leading dimension of the array Z.  LDZ >= max(1,N).
73*
74*  WORK    (workspace) REAL array, dimension (5*N)
75*
76*  IWORK   (workspace) INTEGER array, dimension (N)
77*
78*  IFAIL   (output) INTEGER array, dimension (M)
79*          On normal exit, all elements of IFAIL are zero.
80*          If one or more eigenvectors fail to converge after
81*          MAXITS iterations, then their indices are stored in
82*          array IFAIL.
83*
84*  INFO    (output) INTEGER
85*          = 0: successful exit.
86*          < 0: if INFO = -i, the i-th argument had an illegal value
87*          > 0: if INFO = i, then i eigenvectors failed to converge
88*               in MAXITS iterations.  Their indices are stored in
89*               array IFAIL.
90*
91*  Internal Parameters
92*  ===================
93*
94*  MAXITS  INTEGER, default = 5
95*          The maximum number of iterations performed.
96*
97*  EXTRA   INTEGER, default = 2
98*          The number of iterations performed after norm growth
99*          criterion is satisfied, should be at least 1.
100*
101*  =====================================================================
102*
103*     .. Parameters ..
104      REAL               ZERO, ONE, TEN, ODM3, ODM1
105      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+1,
106     $                   ODM3 = 1.0E-3, ODM1 = 1.0E-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      REAL               CTR, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
115     $                   SCL, SEP, STPCRT, TOL, XJ, XJM
116*     ..
117*     .. Local Arrays ..
118      INTEGER            ISEED( 4 )
119*     ..
120*     .. External Functions ..
121      INTEGER            ISAMAX
122      REAL               SASUM, SDOT, SLAMCH, SNRM2
123      EXTERNAL           ISAMAX, SASUM, SDOT, SLAMCH, SNRM2
124*     ..
125*     .. External Subroutines ..
126      EXTERNAL           SAXPY, SCOPY, SLAGTF, SLAGTS, SLARNV, SSCAL,
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( 'SSTEIN', -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 = SLAMCH( 'Precision' )
179*
180*     Initialize seed for random number generator SLARNV.
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         STPCRT = 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 SLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) )
260*
261*           Copy the matrix T so it won't be destroyed in factorization.
262*
263            CALL SCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 )
264            CALL SCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 )
265            CALL SCOPY( 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 SLAGTF( 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     $            SASUM( BLKSIZ, WORK( INDRV1+1 ), 1 )
286            CALL SSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
287*
288*           Solve the system LU = Pb.
289*
290            CALL SLAGTS( -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                  CTR = -SDOT( BLKSIZ, WORK( INDRV1+1 ), 1, Z( B1, I ),
304     $                  1 )
305                  CALL SAXPY( BLKSIZ, CTR, 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 = ISAMAX( 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.STPCRT )
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 / SNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 )
338            JMAX = ISAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
339            IF( WORK( INDRV1+JMAX ).LT.ZERO )
340     $         SCL = -SCL
341            CALL SSCAL( 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 SSTEIN
361*
362      END
363