1*> \brief \b SSTEDC
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SSTEDC + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sstedc.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sstedc.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sstedc.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
22*                          LIWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          COMPZ
26*       INTEGER            INFO, LDZ, LIWORK, LWORK, N
27*       ..
28*       .. Array Arguments ..
29*       INTEGER            IWORK( * )
30*       REAL               D( * ), E( * ), WORK( * ), Z( LDZ, * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> SSTEDC computes all eigenvalues and, optionally, eigenvectors of a
40*> symmetric tridiagonal matrix using the divide and conquer method.
41*> The eigenvectors of a full or band real symmetric matrix can also be
42*> found if SSYTRD or SSPTRD or SSBTRD has been used to reduce this
43*> matrix to tridiagonal form.
44*>
45*> This code makes very mild assumptions about floating point
46*> arithmetic. It will work on machines with a guard digit in
47*> add/subtract, or on those binary machines without guard digits
48*> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
49*> It could conceivably fail on hexadecimal or decimal machines
50*> without guard digits, but we know of none.  See SLAED3 for details.
51*> \endverbatim
52*
53*  Arguments:
54*  ==========
55*
56*> \param[in] COMPZ
57*> \verbatim
58*>          COMPZ is CHARACTER*1
59*>          = 'N':  Compute eigenvalues only.
60*>          = 'I':  Compute eigenvectors of tridiagonal matrix also.
61*>          = 'V':  Compute eigenvectors of original dense symmetric
62*>                  matrix also.  On entry, Z contains the orthogonal
63*>                  matrix used to reduce the original matrix to
64*>                  tridiagonal form.
65*> \endverbatim
66*>
67*> \param[in] N
68*> \verbatim
69*>          N is INTEGER
70*>          The dimension of the symmetric tridiagonal matrix.  N >= 0.
71*> \endverbatim
72*>
73*> \param[in,out] D
74*> \verbatim
75*>          D is REAL array, dimension (N)
76*>          On entry, the diagonal elements of the tridiagonal matrix.
77*>          On exit, if INFO = 0, the eigenvalues in ascending order.
78*> \endverbatim
79*>
80*> \param[in,out] E
81*> \verbatim
82*>          E is REAL array, dimension (N-1)
83*>          On entry, the subdiagonal elements of the tridiagonal matrix.
84*>          On exit, E has been destroyed.
85*> \endverbatim
86*>
87*> \param[in,out] Z
88*> \verbatim
89*>          Z is REAL array, dimension (LDZ,N)
90*>          On entry, if COMPZ = 'V', then Z contains the orthogonal
91*>          matrix used in the reduction to tridiagonal form.
92*>          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
93*>          orthonormal eigenvectors of the original symmetric matrix,
94*>          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
95*>          of the symmetric tridiagonal matrix.
96*>          If  COMPZ = 'N', then Z is not referenced.
97*> \endverbatim
98*>
99*> \param[in] LDZ
100*> \verbatim
101*>          LDZ is INTEGER
102*>          The leading dimension of the array Z.  LDZ >= 1.
103*>          If eigenvectors are desired, then LDZ >= max(1,N).
104*> \endverbatim
105*>
106*> \param[out] WORK
107*> \verbatim
108*>          WORK is REAL array, dimension (MAX(1,LWORK))
109*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
110*> \endverbatim
111*>
112*> \param[in] LWORK
113*> \verbatim
114*>          LWORK is INTEGER
115*>          The dimension of the array WORK.
116*>          If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
117*>          If COMPZ = 'V' and N > 1 then LWORK must be at least
118*>                         ( 1 + 3*N + 2*N*lg N + 4*N**2 ),
119*>                         where lg( N ) = smallest integer k such
120*>                         that 2**k >= N.
121*>          If COMPZ = 'I' and N > 1 then LWORK must be at least
122*>                         ( 1 + 4*N + N**2 ).
123*>          Note that for COMPZ = 'I' or 'V', then if N is less than or
124*>          equal to the minimum divide size, usually 25, then LWORK need
125*>          only be max(1,2*(N-1)).
126*>
127*>          If LWORK = -1, then a workspace query is assumed; the routine
128*>          only calculates the optimal size of the WORK array, returns
129*>          this value as the first entry of the WORK array, and no error
130*>          message related to LWORK is issued by XERBLA.
131*> \endverbatim
132*>
133*> \param[out] IWORK
134*> \verbatim
135*>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
136*>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
137*> \endverbatim
138*>
139*> \param[in] LIWORK
140*> \verbatim
141*>          LIWORK is INTEGER
142*>          The dimension of the array IWORK.
143*>          If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
144*>          If COMPZ = 'V' and N > 1 then LIWORK must be at least
145*>                         ( 6 + 6*N + 5*N*lg N ).
146*>          If COMPZ = 'I' and N > 1 then LIWORK must be at least
147*>                         ( 3 + 5*N ).
148*>          Note that for COMPZ = 'I' or 'V', then if N is less than or
149*>          equal to the minimum divide size, usually 25, then LIWORK
150*>          need only be 1.
151*>
152*>          If LIWORK = -1, then a workspace query is assumed; the
153*>          routine only calculates the optimal size of the IWORK array,
154*>          returns this value as the first entry of the IWORK array, and
155*>          no error message related to LIWORK is issued by XERBLA.
156*> \endverbatim
157*>
158*> \param[out] INFO
159*> \verbatim
160*>          INFO is INTEGER
161*>          = 0:  successful exit.
162*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
163*>          > 0:  The algorithm failed to compute an eigenvalue while
164*>                working on the submatrix lying in rows and columns
165*>                INFO/(N+1) through mod(INFO,N+1).
166*> \endverbatim
167*
168*  Authors:
169*  ========
170*
171*> \author Univ. of Tennessee
172*> \author Univ. of California Berkeley
173*> \author Univ. of Colorado Denver
174*> \author NAG Ltd.
175*
176*> \date November 2015
177*
178*> \ingroup auxOTHERcomputational
179*
180*> \par Contributors:
181*  ==================
182*>
183*> Jeff Rutter, Computer Science Division, University of California
184*> at Berkeley, USA \n
185*>  Modified by Francoise Tisseur, University of Tennessee
186*>
187*  =====================================================================
188      SUBROUTINE SSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
189     $                   LIWORK, INFO )
190*
191*  -- LAPACK computational routine (version 3.6.0) --
192*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
193*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
194*     November 2015
195*
196*     .. Scalar Arguments ..
197      CHARACTER          COMPZ
198      INTEGER            INFO, LDZ, LIWORK, LWORK, N
199*     ..
200*     .. Array Arguments ..
201      INTEGER            IWORK( * )
202      REAL               D( * ), E( * ), WORK( * ), Z( LDZ, * )
203*     ..
204*
205*  =====================================================================
206*
207*     .. Parameters ..
208      REAL               ZERO, ONE, TWO
209      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
210*     ..
211*     .. Local Scalars ..
212      LOGICAL            LQUERY
213      INTEGER            FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
214     $                   LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
215      REAL               EPS, ORGNRM, P, TINY
216*     ..
217*     .. External Functions ..
218      LOGICAL            LSAME
219      INTEGER            ILAENV
220      REAL               SLAMCH, SLANST
221      EXTERNAL           ILAENV, LSAME, SLAMCH, SLANST
222*     ..
223*     .. External Subroutines ..
224      EXTERNAL           SGEMM, SLACPY, SLAED0, SLASCL, SLASET, SLASRT,
225     $                   SSTEQR, SSTERF, SSWAP, XERBLA
226*     ..
227*     .. Intrinsic Functions ..
228      INTRINSIC          ABS, INT, LOG, MAX, MOD, REAL, SQRT
229*     ..
230*     .. Executable Statements ..
231*
232*     Test the input parameters.
233*
234      INFO = 0
235      LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
236*
237      IF( LSAME( COMPZ, 'N' ) ) THEN
238         ICOMPZ = 0
239      ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
240         ICOMPZ = 1
241      ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
242         ICOMPZ = 2
243      ELSE
244         ICOMPZ = -1
245      END IF
246      IF( ICOMPZ.LT.0 ) THEN
247         INFO = -1
248      ELSE IF( N.LT.0 ) THEN
249         INFO = -2
250      ELSE IF( ( LDZ.LT.1 ) .OR.
251     $         ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
252         INFO = -6
253      END IF
254*
255      IF( INFO.EQ.0 ) THEN
256*
257*        Compute the workspace requirements
258*
259         SMLSIZ = ILAENV( 9, 'SSTEDC', ' ', 0, 0, 0, 0 )
260         IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN
261            LIWMIN = 1
262            LWMIN = 1
263         ELSE IF( N.LE.SMLSIZ ) THEN
264            LIWMIN = 1
265            LWMIN = 2*( N - 1 )
266         ELSE
267            LGN = INT( LOG( REAL( N ) )/LOG( TWO ) )
268            IF( 2**LGN.LT.N )
269     $         LGN = LGN + 1
270            IF( 2**LGN.LT.N )
271     $         LGN = LGN + 1
272            IF( ICOMPZ.EQ.1 ) THEN
273               LWMIN = 1 + 3*N + 2*N*LGN + 4*N**2
274               LIWMIN = 6 + 6*N + 5*N*LGN
275            ELSE IF( ICOMPZ.EQ.2 ) THEN
276               LWMIN = 1 + 4*N + N**2
277               LIWMIN = 3 + 5*N
278            END IF
279         END IF
280         WORK( 1 ) = LWMIN
281         IWORK( 1 ) = LIWMIN
282*
283         IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN
284            INFO = -8
285         ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN
286            INFO = -10
287         END IF
288      END IF
289*
290      IF( INFO.NE.0 ) THEN
291         CALL XERBLA( 'SSTEDC', -INFO )
292         RETURN
293      ELSE IF (LQUERY) THEN
294         RETURN
295      END IF
296*
297*     Quick return if possible
298*
299      IF( N.EQ.0 )
300     $   RETURN
301      IF( N.EQ.1 ) THEN
302         IF( ICOMPZ.NE.0 )
303     $      Z( 1, 1 ) = ONE
304         RETURN
305      END IF
306*
307*     If the following conditional clause is removed, then the routine
308*     will use the Divide and Conquer routine to compute only the
309*     eigenvalues, which requires (3N + 3N**2) real workspace and
310*     (2 + 5N + 2N lg(N)) integer workspace.
311*     Since on many architectures SSTERF is much faster than any other
312*     algorithm for finding eigenvalues only, it is used here
313*     as the default. If the conditional clause is removed, then
314*     information on the size of workspace needs to be changed.
315*
316*     If COMPZ = 'N', use SSTERF to compute the eigenvalues.
317*
318      IF( ICOMPZ.EQ.0 ) THEN
319         CALL SSTERF( N, D, E, INFO )
320         GO TO 50
321      END IF
322*
323*     If N is smaller than the minimum divide size (SMLSIZ+1), then
324*     solve the problem with another solver.
325*
326      IF( N.LE.SMLSIZ ) THEN
327*
328         CALL SSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
329*
330      ELSE
331*
332*        If COMPZ = 'V', the Z matrix must be stored elsewhere for later
333*        use.
334*
335         IF( ICOMPZ.EQ.1 ) THEN
336            STOREZ = 1 + N*N
337         ELSE
338            STOREZ = 1
339         END IF
340*
341         IF( ICOMPZ.EQ.2 ) THEN
342            CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
343         END IF
344*
345*        Scale.
346*
347         ORGNRM = SLANST( 'M', N, D, E )
348         IF( ORGNRM.EQ.ZERO )
349     $      GO TO 50
350*
351         EPS = SLAMCH( 'Epsilon' )
352*
353         START = 1
354*
355*        while ( START <= N )
356*
357   10    CONTINUE
358         IF( START.LE.N ) THEN
359*
360*           Let FINISH be the position of the next subdiagonal entry
361*           such that E( FINISH ) <= TINY or FINISH = N if no such
362*           subdiagonal exists.  The matrix identified by the elements
363*           between START and FINISH constitutes an independent
364*           sub-problem.
365*
366            FINISH = START
367   20       CONTINUE
368            IF( FINISH.LT.N ) THEN
369               TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
370     $                    SQRT( ABS( D( FINISH+1 ) ) )
371               IF( ABS( E( FINISH ) ).GT.TINY ) THEN
372                  FINISH = FINISH + 1
373                  GO TO 20
374               END IF
375            END IF
376*
377*           (Sub) Problem determined.  Compute its size and solve it.
378*
379            M = FINISH - START + 1
380            IF( M.EQ.1 ) THEN
381               START = FINISH + 1
382               GO TO 10
383            END IF
384            IF( M.GT.SMLSIZ ) THEN
385*
386*              Scale.
387*
388               ORGNRM = SLANST( 'M', M, D( START ), E( START ) )
389               CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
390     $                      INFO )
391               CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
392     $                      M-1, INFO )
393*
394               IF( ICOMPZ.EQ.1 ) THEN
395                  STRTRW = 1
396               ELSE
397                  STRTRW = START
398               END IF
399               CALL SLAED0( ICOMPZ, N, M, D( START ), E( START ),
400     $                      Z( STRTRW, START ), LDZ, WORK( 1 ), N,
401     $                      WORK( STOREZ ), IWORK, INFO )
402               IF( INFO.NE.0 ) THEN
403                  INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
404     $                   MOD( INFO, ( M+1 ) ) + START - 1
405                  GO TO 50
406               END IF
407*
408*              Scale back.
409*
410               CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
411     $                      INFO )
412*
413            ELSE
414               IF( ICOMPZ.EQ.1 ) THEN
415*
416*                 Since QR won't update a Z matrix which is larger than
417*                 the length of D, we must solve the sub-problem in a
418*                 workspace and then multiply back into Z.
419*
420                  CALL SSTEQR( 'I', M, D( START ), E( START ), WORK, M,
421     $                         WORK( M*M+1 ), INFO )
422                  CALL SLACPY( 'A', N, M, Z( 1, START ), LDZ,
423     $                         WORK( STOREZ ), N )
424                  CALL SGEMM( 'N', 'N', N, M, M, ONE,
425     $                        WORK( STOREZ ), N, WORK, M, ZERO,
426     $                        Z( 1, START ), LDZ )
427               ELSE IF( ICOMPZ.EQ.2 ) THEN
428                  CALL SSTEQR( 'I', M, D( START ), E( START ),
429     $                         Z( START, START ), LDZ, WORK, INFO )
430               ELSE
431                  CALL SSTERF( M, D( START ), E( START ), INFO )
432               END IF
433               IF( INFO.NE.0 ) THEN
434                  INFO = START*( N+1 ) + FINISH
435                  GO TO 50
436               END IF
437            END IF
438*
439            START = FINISH + 1
440            GO TO 10
441         END IF
442*
443*        endwhile
444*
445         IF( ICOMPZ.EQ.0 ) THEN
446*
447*          Use Quick Sort
448*
449           CALL SLASRT( 'I', N, D, INFO )
450*
451         ELSE
452*
453*          Use Selection Sort to minimize swaps of eigenvectors
454*
455           DO 40 II = 2, N
456              I = II - 1
457              K = I
458              P = D( I )
459              DO 30 J = II, N
460                 IF( D( J ).LT.P ) THEN
461                    K = J
462                    P = D( J )
463                 END IF
464   30         CONTINUE
465              IF( K.NE.I ) THEN
466                 D( K ) = D( I )
467                 D( I ) = P
468                 CALL SSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
469              END IF
470   40      CONTINUE
471         END IF
472      END IF
473*
474   50 CONTINUE
475      WORK( 1 ) = LWMIN
476      IWORK( 1 ) = LIWMIN
477*
478      RETURN
479*
480*     End of SSTEDC
481*
482      END
483