LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date November 2011
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.4.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 2011
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 * If the problem split any number of times, then the eigenvalues
457 * will not be properly ordered. Here we permute the eigenvalues
458 * (and the associated eigenvectors) into ascending order.
459 *
460  IF( m.NE.n ) THEN
461 *
462 * Use Selection Sort to minimize swaps of eigenvectors
463 *
464  DO 60 ii = 2, n
465  i = ii - 1
466  k = i
467  p = d( i )
468  DO 50 j = ii, n
469  IF( d( j ).LT.p ) THEN
470  k = j
471  p = d( j )
472  END IF
473  50 continue
474  IF( k.NE.i ) THEN
475  d( k ) = d( i )
476  d( i ) = p
477  CALL cswap( n, z( 1, i ), 1, z( 1, k ), 1 )
478  END IF
479  60 continue
480  END IF
481  END IF
482 *
483  70 continue
484  work( 1 ) = lwmin
485  rwork( 1 ) = lrwmin
486  iwork( 1 ) = liwmin
487 *
488  return
489 *
490 * End of CSTEDC
491 *
492  END