1      SUBROUTINE CSTEQR2( COMPZ, N, D, E, Z, LDZ, NR, WORK, INFO )
2*
3*  -- ScaLAPACK routine (version 1.7) --
4*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5*     Courant Institute, Argonne National Lab, and Rice University
6*     November 15, 1997
7*
8*     .. Scalar Arguments ..
9      CHARACTER          COMPZ
10      INTEGER            INFO, LDZ, N, NR
11*     ..
12*     .. Array Arguments ..
13      REAL               D( * ), E( * ), WORK( * )
14      COMPLEX            Z( LDZ, * )
15*     ..
16*
17*  Purpose
18*  =======
19*
20*  CSTEQR2 is a modified version of LAPACK routine CSTEQR.
21*  CSTEQR2 computes all eigenvalues and, optionally, eigenvectors of a
22*  symmetric tridiagonal matrix using the implicit QL or QR method.
23*  CSTEQR2 is modified from CSTEQR to allow each ScaLAPACK process
24*  running CSTEQR2 to perform updates on a distributed matrix Q.
25*  Proper usage of CSTEQR2 can be gleaned from
26*  examination of ScaLAPACK's *  PCHEEV.
27*  CSTEQR2 incorporates changes attributed to Greg Henry.
28*
29*  Arguments
30*  =========
31*
32*  COMPZ   (input) CHARACTER*1
33*          = 'N':  Compute eigenvalues only.
34*          = 'I':  Compute eigenvalues and eigenvectors of the
35*                  tridiagonal matrix.  Z must be initialized to the
36*                  identity matrix by PCLASET or CLASET prior
37*                  to entering this subroutine.
38*
39*  N       (input) INTEGER
40*          The order of the matrix.  N >= 0.
41*
42*  D       (input/output) REAL array, dimension (N)
43*          On entry, the diagonal elements of the tridiagonal matrix.
44*          On exit, if INFO = 0, the eigenvalues in ascending order.
45*
46*  E       (input/output) REAL array, dimension (N-1)
47*          On entry, the (n-1) subdiagonal elements of the tridiagonal
48*          matrix.
49*          On exit, E has been destroyed.
50*
51*  Z       (local input/local output) COMPLEX array, global
52*          dimension (N, N), local dimension (LDZ, NR).
53*          On entry, if  COMPZ = 'V', then Z contains the orthogonal
54*          matrix used in the reduction to tridiagonal form.
55*          On exit, if INFO = 0, then if  COMPZ = 'V', Z contains the
56*          orthonormal eigenvectors of the original symmetric matrix,
57*          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
58*          of the symmetric tridiagonal matrix.
59*          If COMPZ = 'N', then Z is not referenced.
60*
61*  LDZ     (input) INTEGER
62*          The leading dimension of the array Z.  LDZ >= 1, and if
63*          eigenvectors are desired, then  LDZ >= max(1,N).
64*
65*  NR      (input) INTEGER
66*          NR = MAX(1, NUMROC( N, NB, MYPROW, 0, NPROCS ) ).
67*          If COMPZ = 'N', then NR is not referenced.
68*
69*  WORK    (workspace) REAL array, dimension (max(1,2*N-2))
70*          If COMPZ = 'N', then WORK is not referenced.
71*
72*  INFO    (output) INTEGER
73*          = 0:  successful exit
74*          < 0:  if INFO = -i, the i-th argument had an illegal value
75*          > 0:  the algorithm has failed to find all the eigenvalues in
76*                a total of 30*N iterations; if INFO = i, then i
77*                elements of E have not converged to zero; on exit, D
78*                and E contain the elements of a symmetric tridiagonal
79*                matrix which is orthogonally similar to the original
80*                matrix.
81*
82*  =====================================================================
83*
84*     .. Parameters ..
85      REAL               ZERO, ONE, TWO, THREE, HALF
86      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
87     $                   THREE = 3.0E0, HALF = 0.5E0 )
88      COMPLEX            CONE
89      PARAMETER          ( CONE = ( 1.0E0, 1.0E0 ) )
90      INTEGER            MAXIT, NMAXLOOK
91      PARAMETER          ( MAXIT = 30, NMAXLOOK = 15 )
92*     ..
93*     .. Local Scalars ..
94      INTEGER            I, ICOMPZ, II, ILAST, ISCALE, J, JTOT, K, L,
95     $                   L1, LEND, LENDM1, LENDP1, LENDSV, LM1, LSV, M,
96     $                   MM, MM1, NLOOK, NM1, NMAXIT
97      REAL               ANORM, B, C, EPS, EPS2, F, G, GP, OLDEL, OLDGP,
98     $                   OLDRP, P, R, RP, RT1, RT2, S, SAFMAX, SAFMIN,
99     $                   SSFMAX, SSFMIN, TST, TST1
100*     ..
101*     .. External Functions ..
102      LOGICAL            LSAME
103      REAL               SLAMCH, SLANST, SLAPY2
104      EXTERNAL           LSAME, SLAMCH, SLANST, SLAPY2
105*     ..
106*     .. External Subroutines ..
107      EXTERNAL           CLASR, CSWAP, SLAEV2, SLARTG, SLASCL, SSTERF,
108     $                   XERBLA
109*     ..
110*     .. Intrinsic Functions ..
111      INTRINSIC          ABS, MAX, SIGN, SQRT
112*     ..
113*     .. Executable Statements ..
114*
115*     Test the input parameters.
116*
117      ILAST = 0
118      INFO = 0
119*
120      IF( LSAME( COMPZ, 'N' ) ) THEN
121         ICOMPZ = 0
122      ELSEIF( LSAME( COMPZ, 'I' ) ) THEN
123         ICOMPZ = 1
124      ELSE
125         ICOMPZ = -1
126      ENDIF
127      IF( ICOMPZ.LT.0 ) THEN
128         INFO = -1
129      ELSEIF( N.LT.0 ) THEN
130         INFO = -2
131      ELSEIF( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, NR ) ) THEN
132         INFO = -6
133      ENDIF
134      IF( INFO.NE.0 ) THEN
135         CALL XERBLA( 'CSTEQR2', -INFO )
136         RETURN
137      ENDIF
138*
139*     Quick return if possible
140*
141      IF( N.EQ.0 )
142     $   RETURN
143*
144*     If eigenvectors aren't not desired, this is faster
145*
146      IF( ICOMPZ.EQ.0 ) THEN
147         CALL SSTERF( N, D, E, INFO )
148         RETURN
149      ENDIF
150*
151      IF( N.EQ.1 ) THEN
152         Z( 1, 1 ) = CONE
153         RETURN
154      ENDIF
155*
156*     Determine the unit roundoff and over/underflow thresholds.
157*
158      EPS = SLAMCH( 'E' )
159      EPS2 = EPS**2
160      SAFMIN = SLAMCH( 'S' )
161      SAFMAX = ONE / SAFMIN
162      SSFMAX = SQRT( SAFMAX ) / THREE
163      SSFMIN = SQRT( SAFMIN ) / EPS2
164*
165*     Compute the eigenvalues and eigenvectors of the tridiagonal
166*     matrix.
167*
168      NMAXIT = N*MAXIT
169      JTOT = 0
170*
171*     Determine where the matrix splits and choose QL or QR iteration
172*     for each block, according to whether top or bottom diagonal
173*     element is smaller.
174*
175      L1 = 1
176      NM1 = N - 1
177*
178   10 CONTINUE
179      IF( L1.GT.N )
180     $   GOTO 220
181      IF( L1.GT.1 )
182     $   E( L1-1 ) = ZERO
183      IF( L1.LE.NM1 ) THEN
184         DO 20 M = L1, NM1
185            TST = ABS( E( M ) )
186            IF( TST.EQ.ZERO )
187     $         GOTO 30
188            IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
189     $          1 ) ) ) )*EPS ) THEN
190               E( M ) = ZERO
191               GOTO 30
192            ENDIF
193   20    CONTINUE
194      ENDIF
195      M = N
196*
197   30 CONTINUE
198      L = L1
199      LSV = L
200      LEND = M
201      LENDSV = LEND
202      L1 = M + 1
203      IF( LEND.EQ.L )
204     $   GOTO 10
205*
206*     Scale submatrix in rows and columns L to LEND
207*
208      ANORM = SLANST( 'I', LEND-L+1, D( L ), E( L ) )
209      ISCALE = 0
210      IF( ANORM.EQ.ZERO )
211     $   GOTO 10
212      IF( ANORM.GT.SSFMAX ) THEN
213         ISCALE = 1
214         CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
215     $                INFO )
216         CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
217     $                INFO )
218      ELSEIF( ANORM.LT.SSFMIN ) THEN
219         ISCALE = 2
220         CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
221     $                INFO )
222         CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
223     $                INFO )
224      ENDIF
225*
226*     Choose between QL and QR iteration
227*
228      IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
229         LEND = LSV
230         L = LENDSV
231      ENDIF
232*
233      IF( LEND.GT.L ) THEN
234*
235*        QL Iteration
236*
237*        Look for small subdiagonal element.
238*
239   40    CONTINUE
240         IF( L.NE.LEND ) THEN
241            LENDM1 = LEND - 1
242            DO 50 M = L, LENDM1
243               TST = ABS( E( M ) )**2
244               IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+
245     $             SAFMIN )GOTO 60
246   50       CONTINUE
247         ENDIF
248*
249         M = LEND
250*
251   60    CONTINUE
252         IF( M.LT.LEND )
253     $      E( M ) = ZERO
254         P = D( L )
255         IF( M.EQ.L )
256     $      GOTO 110
257*
258*        If remaining matrix is 2-by-2, use SLAE2 or SLAEV2
259*        to compute its eigensystem.
260*
261         IF( M.EQ.L+1 ) THEN
262            CALL SLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
263            WORK( L ) = C
264            WORK( N-1+L ) = S
265            CALL CLASR( 'R', 'V', 'B', NR, 2, WORK( L ), WORK( N-1+L ),
266     $                  Z( 1, L ), LDZ )
267            D( L ) = RT1
268            D( L+1 ) = RT2
269            E( L ) = ZERO
270            L = L + 2
271            IF( L.LE.LEND )
272     $         GOTO 40
273            GOTO 200
274         ENDIF
275*
276         IF( JTOT.EQ.NMAXIT )
277     $      GOTO 200
278         JTOT = JTOT + 1
279*
280*        Form shift.
281*
282         G = ( D( L+1 )-P ) / ( TWO*E( L ) )
283         R = SLAPY2( G, ONE )
284         G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) )
285*
286         IF( ICOMPZ.EQ.0 ) THEN
287*           Do not do a lookahead!
288            GOTO 90
289         ENDIF
290*
291         OLDEL = ABS( E( L ) )
292         GP = G
293         RP = R
294         TST = ABS( E( L ) )**2
295         TST = TST / ( ( EPS2*ABS( D( L ) ) )*ABS( D( L+1 ) )+SAFMIN )
296*
297         NLOOK = 1
298         IF( ( TST.GT.ONE ) .AND. ( NLOOK.LE.NMAXLOOK ) ) THEN
299   70       CONTINUE
300*
301*           This is the lookahead loop, going until we have
302*           convergence or too many steps have been taken.
303*
304            S = ONE
305            C = ONE
306            P = ZERO
307            MM1 = M - 1
308            DO 80 I = MM1, L, -1
309               F = S*E( I )
310               B = C*E( I )
311               CALL SLARTG( GP, F, C, S, RP )
312               GP = D( I+1 ) - P
313               RP = ( D( I )-GP )*S + TWO*C*B
314               P = S*RP
315               IF( I.NE.L )
316     $            GP = C*RP - B
317   80       CONTINUE
318            OLDGP = GP
319            OLDRP = RP
320*           Find GP & RP for the next iteration
321            IF( ABS( C*OLDRP-B ).GT.SAFMIN ) THEN
322               GP = ( ( OLDGP+P )-( D( L )-P ) ) / ( TWO*( C*OLDRP-B ) )
323            ELSE
324*
325*           Goto put in by G. Henry to fix ALPHA problem
326*
327               GOTO 90
328*              GP = ( ( OLDGP+P )-( D( L )-P ) ) /
329*    $              ( TWO*( C*OLDRP-B )+SAFMIN )
330            ENDIF
331            RP = SLAPY2( GP, ONE )
332            GP = D( M ) - ( D( L )-P ) +
333     $           ( ( C*OLDRP-B ) / ( GP+SIGN( RP, GP ) ) )
334            TST1 = TST
335            TST = ABS( C*OLDRP-B )**2
336            TST = TST / ( ( EPS2*ABS( D( L )-P ) )*ABS( OLDGP+P )+
337     $            SAFMIN )
338*           Make sure that we are making progress
339            IF( ABS( C*OLDRP-B ).GT.0.9E0*OLDEL ) THEN
340               IF( ABS( C*OLDRP-B ).GT.OLDEL ) THEN
341                  GP = G
342                  RP = R
343               ENDIF
344               TST = HALF
345            ELSE
346               OLDEL = ABS( C*OLDRP-B )
347            ENDIF
348            NLOOK = NLOOK + 1
349            IF( ( TST.GT.ONE ) .AND. ( NLOOK.LE.NMAXLOOK ) )
350     $         GOTO 70
351         ENDIF
352*
353         IF( ( TST.LE.ONE ) .AND. ( TST.NE.HALF ) .AND.
354     $       ( ABS( P ).LT.EPS*ABS( D( L ) ) ) .AND.
355     $       ( ILAST.EQ.L ) .AND. ( ABS( E( L ) )**2.LE.10000.0E0*
356     $       ( ( EPS2*ABS( D( L ) ) )*ABS( D( L+1 ) )+SAFMIN ) ) ) THEN
357*
358*           Skip the current step: the subdiagonal info is just noise.
359*
360            M = L
361            E( M ) = ZERO
362            P = D( L )
363            JTOT = JTOT - 1
364            GOTO 110
365         ENDIF
366         G = GP
367         R = RP
368*
369*        Lookahead over
370*
371   90    CONTINUE
372*
373         S = ONE
374         C = ONE
375         P = ZERO
376*
377*        Inner loop
378*
379         MM1 = M - 1
380         DO 100 I = MM1, L, -1
381            F = S*E( I )
382            B = C*E( I )
383            CALL SLARTG( G, F, C, S, R )
384            IF( I.NE.M-1 )
385     $         E( I+1 ) = R
386            G = D( I+1 ) - P
387            R = ( D( I )-G )*S + TWO*C*B
388            P = S*R
389            D( I+1 ) = G + P
390            G = C*R - B
391*
392*           If eigenvectors are desired, then save rotations.
393*
394            WORK( I ) = C
395            WORK( N-1+I ) = -S
396*
397  100    CONTINUE
398*
399*        If eigenvectors are desired, then apply saved rotations.
400*
401         MM = M - L + 1
402         CALL CLASR( 'R', 'V', 'B', NR, MM, WORK( L ), WORK( N-1+L ),
403     $               Z( 1, L ), LDZ )
404*
405         D( L ) = D( L ) - P
406         E( L ) = G
407         ILAST = L
408         GOTO 40
409*
410*        Eigenvalue found.
411*
412  110    CONTINUE
413         D( L ) = P
414*
415         L = L + 1
416         IF( L.LE.LEND )
417     $      GOTO 40
418         GOTO 200
419*
420      ELSE
421*
422*        QR Iteration
423*
424*        Look for small superdiagonal element.
425*
426  120    CONTINUE
427         IF( L.NE.LEND ) THEN
428            LENDP1 = LEND + 1
429            DO 130 M = L, LENDP1, -1
430               TST = ABS( E( M-1 ) )**2
431               IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+
432     $             SAFMIN )GOTO 140
433  130       CONTINUE
434         ENDIF
435*
436         M = LEND
437*
438  140    CONTINUE
439         IF( M.GT.LEND )
440     $      E( M-1 ) = ZERO
441         P = D( L )
442         IF( M.EQ.L )
443     $      GOTO 190
444*
445*        If remaining matrix is 2-by-2, use SLAE2 or SLAEV2
446*        to compute its eigensystem.
447*
448         IF( M.EQ.L-1 ) THEN
449            CALL SLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
450            WORK( M ) = C
451            WORK( N-1+M ) = S
452            CALL CLASR( 'R', 'V', 'F', NR, 2, WORK( M ), WORK( N-1+M ),
453     $                  Z( 1, L-1 ), LDZ )
454            D( L-1 ) = RT1
455            D( L ) = RT2
456            E( L-1 ) = ZERO
457            L = L - 2
458            IF( L.GE.LEND )
459     $         GOTO 120
460            GOTO 200
461         ENDIF
462*
463         IF( JTOT.EQ.NMAXIT )
464     $      GOTO 200
465         JTOT = JTOT + 1
466*
467*        Form shift.
468*
469         G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
470         R = SLAPY2( G, ONE )
471         G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
472*
473         IF( ICOMPZ.EQ.0 ) THEN
474*           Do not do a lookahead!
475            GOTO 170
476         ENDIF
477*
478         OLDEL = ABS( E( L-1 ) )
479         GP = G
480         RP = R
481         TST = ABS( E( L-1 ) )**2
482         TST = TST / ( ( EPS2*ABS( D( L ) ) )*ABS( D( L-1 ) )+SAFMIN )
483         NLOOK = 1
484         IF( ( TST.GT.ONE ) .AND. ( NLOOK.LE.NMAXLOOK ) ) THEN
485  150       CONTINUE
486*
487*           This is the lookahead loop, going until we have
488*           convergence or too many steps have been taken.
489*
490            S = ONE
491            C = ONE
492            P = ZERO
493*
494*        Inner loop
495*
496            LM1 = L - 1
497            DO 160 I = M, LM1
498               F = S*E( I )
499               B = C*E( I )
500               CALL SLARTG( GP, F, C, S, RP )
501               GP = D( I ) - P
502               RP = ( D( I+1 )-GP )*S + TWO*C*B
503               P = S*RP
504               IF( I.LT.LM1 )
505     $            GP = C*RP - B
506  160       CONTINUE
507            OLDGP = GP
508            OLDRP = RP
509*           Find GP & RP for the next iteration
510            IF( ABS( C*OLDRP-B ).GT.SAFMIN ) THEN
511               GP = ( ( OLDGP+P )-( D( L )-P ) ) / ( TWO*( C*OLDRP-B ) )
512            ELSE
513*
514*           Goto put in by G. Henry to fix ALPHA problem
515*
516               GOTO 170
517*              GP = ( ( OLDGP+P )-( D( L )-P ) ) /
518*    $              ( TWO*( C*OLDRP-B )+SAFMIN )
519            ENDIF
520            RP = SLAPY2( GP, ONE )
521            GP = D( M ) - ( D( L )-P ) +
522     $           ( ( C*OLDRP-B ) / ( GP+SIGN( RP, GP ) ) )
523            TST1 = TST
524            TST = ABS( ( C*OLDRP-B ) )**2
525            TST = TST / ( ( EPS2*ABS( D( L )-P ) )*ABS( OLDGP+P )+
526     $            SAFMIN )
527*           Make sure that we are making progress
528            IF( ABS( C*OLDRP-B ).GT.0.9E0*OLDEL ) THEN
529               IF( ABS( C*OLDRP-B ).GT.OLDEL ) THEN
530                  GP = G
531                  RP = R
532               ENDIF
533               TST = HALF
534            ELSE
535               OLDEL = ABS( C*OLDRP-B )
536            ENDIF
537            NLOOK = NLOOK + 1
538            IF( ( TST.GT.ONE ) .AND. ( NLOOK.LE.NMAXLOOK ) )
539     $         GOTO 150
540         ENDIF
541         IF( ( TST.LE.ONE ) .AND. ( TST.NE.HALF ) .AND.
542     $       ( ABS( P ).LT.EPS*ABS( D( L ) ) ) .AND.
543     $       ( ILAST.EQ.L ) .AND. ( ABS( E( L-1 ) )**2.LE.10000.0E0*
544     $       ( ( EPS2*ABS( D( L-1 ) ) )*ABS( D( L ) )+SAFMIN ) ) ) THEN
545*
546*           Skip the current step: the subdiagonal info is just noise.
547*
548            M = L
549            E( M-1 ) = ZERO
550            P = D( L )
551            JTOT = JTOT - 1
552            GOTO 190
553         ENDIF
554*
555         G = GP
556         R = RP
557*
558*        Lookahead over
559*
560  170    CONTINUE
561*
562         S = ONE
563         C = ONE
564         P = ZERO
565         DO 180 I = M, LM1
566            F = S*E( I )
567            B = C*E( I )
568            CALL SLARTG( G, F, C, S, R )
569            IF( I.NE.M )
570     $         E( I-1 ) = R
571            G = D( I ) - P
572            R = ( D( I+1 )-G )*S + TWO*C*B
573            P = S*R
574            D( I ) = G + P
575            G = C*R - B
576*
577*           If eigenvectors are desired, then save rotations.
578*
579            WORK( I ) = C
580            WORK( N-1+I ) = S
581*
582  180    CONTINUE
583*
584*        If eigenvectors are desired, then apply saved rotations.
585*
586         MM = L - M + 1
587         CALL CLASR( 'R', 'V', 'F', NR, MM, WORK( M ), WORK( N-1+M ),
588     $               Z( 1, M ), LDZ )
589*
590         D( L ) = D( L ) - P
591         E( LM1 ) = G
592         ILAST = L
593         GOTO 120
594*
595*        Eigenvalue found.
596*
597  190    CONTINUE
598         D( L ) = P
599*
600         L = L - 1
601         IF( L.GE.LEND )
602     $      GOTO 120
603         GOTO 200
604*
605      ENDIF
606*
607*     Undo scaling if necessary
608*
609  200 CONTINUE
610      IF( ISCALE.EQ.1 ) THEN
611         CALL SLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
612     $                D( LSV ), N, INFO )
613         CALL SLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ),
614     $                N, INFO )
615      ELSEIF( ISCALE.EQ.2 ) THEN
616         CALL SLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
617     $                D( LSV ), N, INFO )
618         CALL SLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ),
619     $                N, INFO )
620      ENDIF
621*
622*     Check for no convergence to an eigenvalue after a total
623*     of N*MAXIT iterations.
624*
625      IF( JTOT.LT.NMAXIT )
626     $   GOTO 10
627      DO 210 I = 1, N - 1
628         IF( E( I ).NE.ZERO )
629     $      INFO = INFO + 1
630  210 CONTINUE
631      GOTO 250
632*
633*     Order eigenvalues and eigenvectors.
634*
635  220 CONTINUE
636*
637*        Use Selection Sort to minimize swaps of eigenvectors
638*
639      DO 240 II = 2, N
640         I = II - 1
641         K = I
642         P = D( I )
643         DO 230 J = II, N
644            IF( D( J ).LT.P ) THEN
645               K = J
646               P = D( J )
647            ENDIF
648  230    CONTINUE
649         IF( K.NE.I ) THEN
650            D( K ) = D( I )
651            D( I ) = P
652            CALL CSWAP( NR, Z( 1, I ), 1, Z( 1, K ), 1 )
653         ENDIF
654  240 CONTINUE
655*
656  250 CONTINUE
657*     WRITE( *, FMT = * )'JTOT', JTOT
658      RETURN
659*
660*     End of SSTEQR2
661*
662      END
663