LAPACK  3.6.1 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
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 *> \date November 2015
202 *
203 *> \ingroup complexOTHERcomputational
204 *
205 *> \par Contributors:
206 * ==================
207 *>
208 *> Jeff Rutter, Computer Science Division, University of California
209 *> at Berkeley, USA
210 *
211 * =====================================================================
212  SUBROUTINE cstedc( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK,
213  \$ lrwork, iwork, liwork, info )
214 *
215 * -- LAPACK computational routine (version 3.6.0) --
216 * -- LAPACK is a software package provided by Univ. of Tennessee, --
217 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
218 * November 2015
219 *
220 * .. Scalar Arguments ..
221  CHARACTER COMPZ
222  INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N
223 * ..
224 * .. Array Arguments ..
225  INTEGER IWORK( * )
226  REAL D( * ), E( * ), RWORK( * )
227  COMPLEX WORK( * ), Z( ldz, * )
228 * ..
229 *
230 * =====================================================================
231 *
232 * .. Parameters ..
233  REAL ZERO, ONE, TWO
234  parameter ( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
235 * ..
236 * .. Local Scalars ..
237  LOGICAL LQUERY
238  INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN, LL,
239  \$ lrwmin, lwmin, m, smlsiz, start
240  REAL EPS, ORGNRM, P, TINY
241 * ..
242 * .. External Functions ..
243  LOGICAL LSAME
244  INTEGER ILAENV
245  REAL SLAMCH, SLANST
246  EXTERNAL ilaenv, lsame, slamch, slanst
247 * ..
248 * .. External Subroutines ..
249  EXTERNAL xerbla, clacpy, clacrm, claed0, csteqr, cswap,
251 * ..
252 * .. Intrinsic Functions ..
253  INTRINSIC abs, int, log, max, mod, REAL, SQRT
254 * ..
255 * .. Executable Statements ..
256 *
257 * Test the input parameters.
258 *
259  info = 0
260  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
261 *
262  IF( lsame( compz, 'N' ) ) THEN
263  icompz = 0
264  ELSE IF( lsame( compz, 'V' ) ) THEN
265  icompz = 1
266  ELSE IF( lsame( compz, 'I' ) ) THEN
267  icompz = 2
268  ELSE
269  icompz = -1
270  END IF
271  IF( icompz.LT.0 ) THEN
272  info = -1
273  ELSE IF( n.LT.0 ) THEN
274  info = -2
275  ELSE IF( ( ldz.LT.1 ) .OR.
276  \$ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) ) THEN
277  info = -6
278  END IF
279 *
280  IF( info.EQ.0 ) THEN
281 *
282 * Compute the workspace requirements
283 *
284  smlsiz = ilaenv( 9, 'CSTEDC', ' ', 0, 0, 0, 0 )
285  IF( n.LE.1 .OR. icompz.EQ.0 ) THEN
286  lwmin = 1
287  liwmin = 1
288  lrwmin = 1
289  ELSE IF( n.LE.smlsiz ) THEN
290  lwmin = 1
291  liwmin = 1
292  lrwmin = 2*( n - 1 )
293  ELSE IF( icompz.EQ.1 ) THEN
294  lgn = int( log( REAL( N ) ) / log( TWO ) )
295  IF( 2**lgn.LT.n )
296  \$ lgn = lgn + 1
297  IF( 2**lgn.LT.n )
298  \$ lgn = lgn + 1
299  lwmin = n*n
300  lrwmin = 1 + 3*n + 2*n*lgn + 4*n**2
301  liwmin = 6 + 6*n + 5*n*lgn
302  ELSE IF( icompz.EQ.2 ) THEN
303  lwmin = 1
304  lrwmin = 1 + 4*n + 2*n**2
305  liwmin = 3 + 5*n
306  END IF
307  work( 1 ) = lwmin
308  rwork( 1 ) = lrwmin
309  iwork( 1 ) = liwmin
310 *
311  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
312  info = -8
313  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
314  info = -10
315  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
316  info = -12
317  END IF
318  END IF
319 *
320  IF( info.NE.0 ) THEN
321  CALL xerbla( 'CSTEDC', -info )
322  RETURN
323  ELSE IF( lquery ) THEN
324  RETURN
325  END IF
326 *
327 * Quick return if possible
328 *
329  IF( n.EQ.0 )
330  \$ RETURN
331  IF( n.EQ.1 ) THEN
332  IF( icompz.NE.0 )
333  \$ z( 1, 1 ) = one
334  RETURN
335  END IF
336 *
337 * If the following conditional clause is removed, then the routine
338 * will use the Divide and Conquer routine to compute only the
339 * eigenvalues, which requires (3N + 3N**2) real workspace and
340 * (2 + 5N + 2N lg(N)) integer workspace.
341 * Since on many architectures SSTERF is much faster than any other
342 * algorithm for finding eigenvalues only, it is used here
343 * as the default. If the conditional clause is removed, then
344 * information on the size of workspace needs to be changed.
345 *
346 * If COMPZ = 'N', use SSTERF to compute the eigenvalues.
347 *
348  IF( icompz.EQ.0 ) THEN
349  CALL ssterf( n, d, e, info )
350  GO TO 70
351  END IF
352 *
353 * If N is smaller than the minimum divide size (SMLSIZ+1), then
354 * solve the problem with another solver.
355 *
356  IF( n.LE.smlsiz ) THEN
357 *
358  CALL csteqr( compz, n, d, e, z, ldz, rwork, info )
359 *
360  ELSE
361 *
362 * If COMPZ = 'I', we simply call SSTEDC instead.
363 *
364  IF( icompz.EQ.2 ) THEN
365  CALL slaset( 'Full', n, n, zero, one, rwork, n )
366  ll = n*n + 1
367  CALL sstedc( 'I', n, d, e, rwork, n,
368  \$ rwork( ll ), lrwork-ll+1, iwork, liwork, info )
369  DO 20 j = 1, n
370  DO 10 i = 1, n
371  z( i, j ) = rwork( ( j-1 )*n+i )
372  10 CONTINUE
373  20 CONTINUE
374  GO TO 70
375  END IF
376 *
377 * From now on, only option left to be handled is COMPZ = 'V',
378 * i.e. ICOMPZ = 1.
379 *
380 * Scale.
381 *
382  orgnrm = slanst( 'M', n, d, e )
383  IF( orgnrm.EQ.zero )
384  \$ GO TO 70
385 *
386  eps = slamch( 'Epsilon' )
387 *
388  start = 1
389 *
390 * while ( START <= N )
391 *
392  30 CONTINUE
393  IF( start.LE.n ) THEN
394 *
395 * Let FINISH be the position of the next subdiagonal entry
396 * such that E( FINISH ) <= TINY or FINISH = N if no such
397 * subdiagonal exists. The matrix identified by the elements
398 * between START and FINISH constitutes an independent
399 * sub-problem.
400 *
401  finish = start
402  40 CONTINUE
403  IF( finish.LT.n ) THEN
404  tiny = eps*sqrt( abs( d( finish ) ) )*
405  \$ sqrt( abs( d( finish+1 ) ) )
406  IF( abs( e( finish ) ).GT.tiny ) THEN
407  finish = finish + 1
408  GO TO 40
409  END IF
410  END IF
411 *
412 * (Sub) Problem determined. Compute its size and solve it.
413 *
414  m = finish - start + 1
415  IF( m.GT.smlsiz ) THEN
416 *
417 * Scale.
418 *
419  orgnrm = slanst( 'M', m, d( start ), e( start ) )
420  CALL slascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
421  \$ info )
422  CALL slascl( 'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
423  \$ m-1, info )
424 *
425  CALL claed0( n, m, d( start ), e( start ), z( 1, start ),
426  \$ ldz, work, n, rwork, iwork, info )
427  IF( info.GT.0 ) THEN
428  info = ( info / ( m+1 )+start-1 )*( n+1 ) +
429  \$ mod( info, ( m+1 ) ) + start - 1
430  GO TO 70
431  END IF
432 *
433 * Scale back.
434 *
435  CALL slascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
436  \$ info )
437 *
438  ELSE
439  CALL ssteqr( 'I', m, d( start ), e( start ), rwork, m,
440  \$ rwork( m*m+1 ), info )
441  CALL clacrm( n, m, z( 1, start ), ldz, rwork, m, work, n,
442  \$ rwork( m*m+1 ) )
443  CALL clacpy( 'A', n, m, work, n, z( 1, start ), ldz )
444  IF( info.GT.0 ) THEN
445  info = start*( n+1 ) + finish
446  GO TO 70
447  END IF
448  END IF
449 *
450  start = finish + 1
451  GO TO 30
452  END IF
453 *
454 * endwhile
455 *
456 *
457 * Use Selection Sort to minimize swaps of eigenvectors
458 *
459  DO 60 ii = 2, n
460  i = ii - 1
461  k = i
462  p = d( i )
463  DO 50 j = ii, n
464  IF( d( j ).LT.p ) THEN
465  k = j
466  p = d( j )
467  END IF
468  50 CONTINUE
469  IF( k.NE.i ) THEN
470  d( k ) = d( i )
471  d( i ) = p
472  CALL cswap( n, z( 1, i ), 1, z( 1, k ), 1 )
473  END IF
474  60 CONTINUE
475  END IF
476 *
477  70 CONTINUE
478  work( 1 ) = lwmin
479  rwork( 1 ) = lrwmin
480  iwork( 1 ) = liwmin
481 *
482  RETURN
483 *
484 * End of CSTEDC
485 *
486  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:145
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine csteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
CSTEQR
Definition: csteqr.f:134
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:112
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:133
subroutine claed0(QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK, IWORK, INFO)
CLAED0 used by sstedc. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
Definition: claed0.f:147
subroutine cstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
CSTEDC
Definition: cstedc.f:214
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
subroutine clacrm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
CLACRM multiplies a complex matrix by a square real matrix.
Definition: clacrm.f:116
subroutine sstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEDC
Definition: sstedc.f:190