LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
cstedc.f
Go to the documentation of this file.
1 *> \brief \b CSTEDC
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CSTEDC + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cstedc.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cstedc.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cstedc.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK,
22 * LRWORK, IWORK, LIWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER COMPZ
26 * INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IWORK( * )
30 * REAL D( * ), E( * ), RWORK( * )
31 * COMPLEX WORK( * ), Z( LDZ, * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> CSTEDC computes all eigenvalues and, optionally, eigenvectors of a
41 *> symmetric tridiagonal matrix using the divide and conquer method.
42 *> The eigenvectors of a full or band complex Hermitian matrix can also
43 *> be found if CHETRD or CHPTRD or CHBTRD has been used to reduce this
44 *> matrix to tridiagonal form.
45 *>
46 *> This code makes very mild assumptions about floating point
47 *> arithmetic. It will work on machines with a guard digit in
48 *> add/subtract, or on those binary machines without guard digits
49 *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
50 *> It could conceivably fail on hexadecimal or decimal machines
51 *> without guard digits, but we know of none. See SLAED3 for details.
52 *> \endverbatim
53 *
54 * Arguments:
55 * ==========
56 *
57 *> \param[in] COMPZ
58 *> \verbatim
59 *> COMPZ is CHARACTER*1
60 *> = 'N': Compute eigenvalues only.
61 *> = 'I': Compute eigenvectors of tridiagonal matrix also.
62 *> = 'V': Compute eigenvectors of original Hermitian matrix
63 *> also. On entry, Z contains the unitary matrix used
64 *> to reduce the original matrix to 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 COMPLEX array, dimension (LDZ,N)
90 *> On entry, if COMPZ = 'V', then Z contains the unitary
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 Hermitian 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 COMPLEX 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 'I', or N <= 1, LWORK must be at least 1.
117 *> If COMPZ = 'V' and N > 1, LWORK must be at least N*N.
118 *> Note that for COMPZ = 'V', then if N is less than or
119 *> equal to the minimum divide size, usually 25, then LWORK need
120 *> only be 1.
121 *>
122 *> If LWORK = -1, then a workspace query is assumed; the routine
123 *> only calculates the optimal sizes of the WORK, RWORK and
124 *> IWORK arrays, returns these values as the first entries of
125 *> the WORK, RWORK and IWORK arrays, and no error message
126 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
127 *> \endverbatim
128 *>
129 *> \param[out] RWORK
130 *> \verbatim
131 *> RWORK is REAL array, dimension (MAX(1,LRWORK))
132 *> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
133 *> \endverbatim
134 *>
135 *> \param[in] LRWORK
136 *> \verbatim
137 *> LRWORK is INTEGER
138 *> The dimension of the array RWORK.
139 *> If COMPZ = 'N' or N <= 1, LRWORK must be at least 1.
140 *> If COMPZ = 'V' and N > 1, LRWORK must be at least
141 *> 1 + 3*N + 2*N*lg N + 4*N**2 ,
142 *> where lg( N ) = smallest integer k such
143 *> that 2**k >= N.
144 *> If COMPZ = 'I' and N > 1, LRWORK must be at least
145 *> 1 + 4*N + 2*N**2 .
146 *> Note that for COMPZ = 'I' or 'V', then if N is less than or
147 *> equal to the minimum divide size, usually 25, then LRWORK
148 *> need only be max(1,2*(N-1)).
149 *>
150 *> If LRWORK = -1, then a workspace query is assumed; the
151 *> routine only calculates the optimal sizes of the WORK, RWORK
152 *> and IWORK arrays, returns these values as the first entries
153 *> of the WORK, RWORK and IWORK arrays, and no error message
154 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
155 *> \endverbatim
156 *>
157 *> \param[out] IWORK
158 *> \verbatim
159 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
160 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
161 *> \endverbatim
162 *>
163 *> \param[in] LIWORK
164 *> \verbatim
165 *> LIWORK is INTEGER
166 *> The dimension of the array IWORK.
167 *> If COMPZ = 'N' or N <= 1, LIWORK must be at least 1.
168 *> If COMPZ = 'V' or N > 1, LIWORK must be at least
169 *> 6 + 6*N + 5*N*lg N.
170 *> If COMPZ = 'I' or N > 1, LIWORK must be at least
171 *> 3 + 5*N .
172 *> Note that for COMPZ = 'I' or 'V', then if N is less than or
173 *> equal to the minimum divide size, usually 25, then LIWORK
174 *> need only be 1.
175 *>
176 *> If LIWORK = -1, then a workspace query is assumed; the
177 *> routine only calculates the optimal sizes of the WORK, RWORK
178 *> and IWORK arrays, returns these values as the first entries
179 *> of the WORK, RWORK and IWORK arrays, and no error message
180 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
181 *> \endverbatim
182 *>
183 *> \param[out] INFO
184 *> \verbatim
185 *> INFO is INTEGER
186 *> = 0: successful exit.
187 *> < 0: if INFO = -i, the i-th argument had an illegal value.
188 *> > 0: The algorithm failed to compute an eigenvalue while
189 *> working on the submatrix lying in rows and columns
190 *> INFO/(N+1) through mod(INFO,N+1).
191 *> \endverbatim
192 *
193 * Authors:
194 * ========
195 *
196 *> \author Univ. of Tennessee
197 *> \author Univ. of California Berkeley
198 *> \author Univ. of Colorado Denver
199 *> \author NAG Ltd.
200 *
201 *> \ingroup complexOTHERcomputational
202 *
203 *> \par Contributors:
204 * ==================
205 *>
206 *> Jeff Rutter, Computer Science Division, University of California
207 *> at Berkeley, USA
208 *
209 * =====================================================================
210  SUBROUTINE cstedc( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK,
211  $ LRWORK, IWORK, LIWORK, INFO )
212 *
213 * -- LAPACK computational routine --
214 * -- LAPACK is a software package provided by Univ. of Tennessee, --
215 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
216 *
217 * .. Scalar Arguments ..
218  CHARACTER COMPZ
219  INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N
220 * ..
221 * .. Array Arguments ..
222  INTEGER IWORK( * )
223  REAL D( * ), E( * ), RWORK( * )
224  COMPLEX WORK( * ), Z( LDZ, * )
225 * ..
226 *
227 * =====================================================================
228 *
229 * .. Parameters ..
230  REAL ZERO, ONE, TWO
231  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
232 * ..
233 * .. Local Scalars ..
234  LOGICAL LQUERY
235  INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN, LL,
236  $ lrwmin, lwmin, m, smlsiz, start
237  REAL EPS, ORGNRM, P, TINY
238 * ..
239 * .. External Functions ..
240  LOGICAL LSAME
241  INTEGER ILAENV
242  REAL SLAMCH, SLANST
243  EXTERNAL ilaenv, lsame, slamch, slanst
244 * ..
245 * .. External Subroutines ..
246  EXTERNAL xerbla, clacpy, clacrm, claed0, csteqr, cswap,
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC abs, int, log, max, mod, real, sqrt
251 * ..
252 * .. Executable Statements ..
253 *
254 * Test the input parameters.
255 *
256  info = 0
257  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
258 *
259  IF( lsame( compz, 'N' ) ) THEN
260  icompz = 0
261  ELSE IF( lsame( compz, 'V' ) ) THEN
262  icompz = 1
263  ELSE IF( lsame( compz, 'I' ) ) THEN
264  icompz = 2
265  ELSE
266  icompz = -1
267  END IF
268  IF( icompz.LT.0 ) THEN
269  info = -1
270  ELSE IF( n.LT.0 ) THEN
271  info = -2
272  ELSE IF( ( ldz.LT.1 ) .OR.
273  $ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) ) THEN
274  info = -6
275  END IF
276 *
277  IF( info.EQ.0 ) THEN
278 *
279 * Compute the workspace requirements
280 *
281  smlsiz = ilaenv( 9, 'CSTEDC', ' ', 0, 0, 0, 0 )
282  IF( n.LE.1 .OR. icompz.EQ.0 ) THEN
283  lwmin = 1
284  liwmin = 1
285  lrwmin = 1
286  ELSE IF( n.LE.smlsiz ) THEN
287  lwmin = 1
288  liwmin = 1
289  lrwmin = 2*( n - 1 )
290  ELSE IF( icompz.EQ.1 ) THEN
291  lgn = int( log( real( n ) ) / log( two ) )
292  IF( 2**lgn.LT.n )
293  $ lgn = lgn + 1
294  IF( 2**lgn.LT.n )
295  $ lgn = lgn + 1
296  lwmin = n*n
297  lrwmin = 1 + 3*n + 2*n*lgn + 4*n**2
298  liwmin = 6 + 6*n + 5*n*lgn
299  ELSE IF( icompz.EQ.2 ) THEN
300  lwmin = 1
301  lrwmin = 1 + 4*n + 2*n**2
302  liwmin = 3 + 5*n
303  END IF
304  work( 1 ) = lwmin
305  rwork( 1 ) = lrwmin
306  iwork( 1 ) = liwmin
307 *
308  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
309  info = -8
310  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
311  info = -10
312  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
313  info = -12
314  END IF
315  END IF
316 *
317  IF( info.NE.0 ) THEN
318  CALL xerbla( 'CSTEDC', -info )
319  RETURN
320  ELSE IF( lquery ) THEN
321  RETURN
322  END IF
323 *
324 * Quick return if possible
325 *
326  IF( n.EQ.0 )
327  $ RETURN
328  IF( n.EQ.1 ) THEN
329  IF( icompz.NE.0 )
330  $ z( 1, 1 ) = one
331  RETURN
332  END IF
333 *
334 * If the following conditional clause is removed, then the routine
335 * will use the Divide and Conquer routine to compute only the
336 * eigenvalues, which requires (3N + 3N**2) real workspace and
337 * (2 + 5N + 2N lg(N)) integer workspace.
338 * Since on many architectures SSTERF is much faster than any other
339 * algorithm for finding eigenvalues only, it is used here
340 * as the default. If the conditional clause is removed, then
341 * information on the size of workspace needs to be changed.
342 *
343 * If COMPZ = 'N', use SSTERF to compute the eigenvalues.
344 *
345  IF( icompz.EQ.0 ) THEN
346  CALL ssterf( n, d, e, info )
347  GO TO 70
348  END IF
349 *
350 * If N is smaller than the minimum divide size (SMLSIZ+1), then
351 * solve the problem with another solver.
352 *
353  IF( n.LE.smlsiz ) THEN
354 *
355  CALL csteqr( compz, n, d, e, z, ldz, rwork, info )
356 *
357  ELSE
358 *
359 * If COMPZ = 'I', we simply call SSTEDC instead.
360 *
361  IF( icompz.EQ.2 ) THEN
362  CALL slaset( 'Full', n, n, zero, one, rwork, n )
363  ll = n*n + 1
364  CALL sstedc( 'I', n, d, e, rwork, n,
365  $ rwork( ll ), lrwork-ll+1, iwork, liwork, info )
366  DO 20 j = 1, n
367  DO 10 i = 1, n
368  z( i, j ) = rwork( ( j-1 )*n+i )
369  10 CONTINUE
370  20 CONTINUE
371  GO TO 70
372  END IF
373 *
374 * From now on, only option left to be handled is COMPZ = 'V',
375 * i.e. ICOMPZ = 1.
376 *
377 * Scale.
378 *
379  orgnrm = slanst( 'M', n, d, e )
380  IF( orgnrm.EQ.zero )
381  $ GO TO 70
382 *
383  eps = slamch( 'Epsilon' )
384 *
385  start = 1
386 *
387 * while ( START <= N )
388 *
389  30 CONTINUE
390  IF( start.LE.n ) THEN
391 *
392 * Let FINISH be the position of the next subdiagonal entry
393 * such that E( FINISH ) <= TINY or FINISH = N if no such
394 * subdiagonal exists. The matrix identified by the elements
395 * between START and FINISH constitutes an independent
396 * sub-problem.
397 *
398  finish = start
399  40 CONTINUE
400  IF( finish.LT.n ) THEN
401  tiny = eps*sqrt( abs( d( finish ) ) )*
402  $ sqrt( abs( d( finish+1 ) ) )
403  IF( abs( e( finish ) ).GT.tiny ) THEN
404  finish = finish + 1
405  GO TO 40
406  END IF
407  END IF
408 *
409 * (Sub) Problem determined. Compute its size and solve it.
410 *
411  m = finish - start + 1
412  IF( m.GT.smlsiz ) THEN
413 *
414 * Scale.
415 *
416  orgnrm = slanst( 'M', m, d( start ), e( start ) )
417  CALL slascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
418  $ info )
419  CALL slascl( 'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
420  $ m-1, info )
421 *
422  CALL claed0( n, m, d( start ), e( start ), z( 1, start ),
423  $ ldz, work, n, rwork, iwork, info )
424  IF( info.GT.0 ) THEN
425  info = ( info / ( m+1 )+start-1 )*( n+1 ) +
426  $ mod( info, ( m+1 ) ) + start - 1
427  GO TO 70
428  END IF
429 *
430 * Scale back.
431 *
432  CALL slascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
433  $ info )
434 *
435  ELSE
436  CALL ssteqr( 'I', m, d( start ), e( start ), rwork, m,
437  $ rwork( m*m+1 ), info )
438  CALL clacrm( n, m, z( 1, start ), ldz, rwork, m, work, n,
439  $ rwork( m*m+1 ) )
440  CALL clacpy( 'A', n, m, work, n, z( 1, start ), ldz )
441  IF( info.GT.0 ) THEN
442  info = start*( n+1 ) + finish
443  GO TO 70
444  END IF
445  END IF
446 *
447  start = finish + 1
448  GO TO 30
449  END IF
450 *
451 * endwhile
452 *
453 *
454 * Use Selection Sort to minimize swaps of eigenvectors
455 *
456  DO 60 ii = 2, n
457  i = ii - 1
458  k = i
459  p = d( i )
460  DO 50 j = ii, n
461  IF( d( j ).LT.p ) THEN
462  k = j
463  p = d( j )
464  END IF
465  50 CONTINUE
466  IF( k.NE.i ) THEN
467  d( k ) = d( i )
468  d( i ) = p
469  CALL cswap( n, z( 1, i ), 1, z( 1, k ), 1 )
470  END IF
471  60 CONTINUE
472  END IF
473 *
474  70 CONTINUE
475  work( 1 ) = lwmin
476  rwork( 1 ) = lrwmin
477  iwork( 1 ) = liwmin
478 *
479  RETURN
480 *
481 * End of CSTEDC
482 *
483  END
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:143
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:131
subroutine sstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEDC
Definition: sstedc.f:188
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine clacrm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
CLACRM multiplies a complex matrix by a square real matrix.
Definition: clacrm.f:114
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
subroutine claed0(QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK, IWORK, INFO)
CLAED0 used by CSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
Definition: claed0.f:145
subroutine csteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
CSTEQR
Definition: csteqr.f:132
subroutine cstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
CSTEDC
Definition: cstedc.f:212