LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zstedc.f
Go to the documentation of this file.
1 *> \brief \b ZSTEDC
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZSTEDC + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstedc.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstedc.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstedc.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZSTEDC( 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 * DOUBLE PRECISION D( * ), E( * ), RWORK( * )
31 * COMPLEX*16 WORK( * ), Z( LDZ, * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> ZSTEDC 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 ZHETRD or ZHPTRD or ZHBTRD 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 DLAED3 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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*16 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 DOUBLE PRECISION array,
132 *> dimension (LRWORK)
133 *> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
134 *> \endverbatim
135 *>
136 *> \param[in] LRWORK
137 *> \verbatim
138 *> LRWORK is INTEGER
139 *> The dimension of the array RWORK.
140 *> If COMPZ = 'N' or N <= 1, LRWORK must be at least 1.
141 *> If COMPZ = 'V' and N > 1, LRWORK must be at least
142 *> 1 + 3*N + 2*N*lg N + 4*N**2 ,
143 *> where lg( N ) = smallest integer k such
144 *> that 2**k >= N.
145 *> If COMPZ = 'I' and N > 1, LRWORK must be at least
146 *> 1 + 4*N + 2*N**2 .
147 *> Note that for COMPZ = 'I' or 'V', then if N is less than or
148 *> equal to the minimum divide size, usually 25, then LRWORK
149 *> need only be max(1,2*(N-1)).
150 *>
151 *> If LRWORK = -1, then a workspace query is assumed; the
152 *> routine only calculates the optimal sizes of the WORK, RWORK
153 *> and IWORK arrays, returns these values as the first entries
154 *> of the WORK, RWORK and IWORK arrays, and no error message
155 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
156 *> \endverbatim
157 *>
158 *> \param[out] IWORK
159 *> \verbatim
160 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
161 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
162 *> \endverbatim
163 *>
164 *> \param[in] LIWORK
165 *> \verbatim
166 *> LIWORK is INTEGER
167 *> The dimension of the array IWORK.
168 *> If COMPZ = 'N' or N <= 1, LIWORK must be at least 1.
169 *> If COMPZ = 'V' or N > 1, LIWORK must be at least
170 *> 6 + 6*N + 5*N*lg N.
171 *> If COMPZ = 'I' or N > 1, LIWORK must be at least
172 *> 3 + 5*N .
173 *> Note that for COMPZ = 'I' or 'V', then if N is less than or
174 *> equal to the minimum divide size, usually 25, then LIWORK
175 *> need only be 1.
176 *>
177 *> If LIWORK = -1, then a workspace query is assumed; the
178 *> routine only calculates the optimal sizes of the WORK, RWORK
179 *> and IWORK arrays, returns these values as the first entries
180 *> of the WORK, RWORK and IWORK arrays, and no error message
181 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
182 *> \endverbatim
183 *>
184 *> \param[out] INFO
185 *> \verbatim
186 *> INFO is INTEGER
187 *> = 0: successful exit.
188 *> < 0: if INFO = -i, the i-th argument had an illegal value.
189 *> > 0: The algorithm failed to compute an eigenvalue while
190 *> working on the submatrix lying in rows and columns
191 *> INFO/(N+1) through mod(INFO,N+1).
192 *> \endverbatim
193 *
194 * Authors:
195 * ========
196 *
197 *> \author Univ. of Tennessee
198 *> \author Univ. of California Berkeley
199 *> \author Univ. of Colorado Denver
200 *> \author NAG Ltd.
201 *
202 *> \date November 2011
203 *
204 *> \ingroup complex16OTHERcomputational
205 *
206 *> \par Contributors:
207 * ==================
208 *>
209 *> Jeff Rutter, Computer Science Division, University of California
210 *> at Berkeley, USA
211 *
212 * =====================================================================
213  SUBROUTINE zstedc( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK,
214  $ lrwork, iwork, liwork, info )
215 *
216 * -- LAPACK computational routine (version 3.4.0) --
217 * -- LAPACK is a software package provided by Univ. of Tennessee, --
218 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
219 * November 2011
220 *
221 * .. Scalar Arguments ..
222  CHARACTER compz
223  INTEGER info, ldz, liwork, lrwork, lwork, n
224 * ..
225 * .. Array Arguments ..
226  INTEGER iwork( * )
227  DOUBLE PRECISION d( * ), e( * ), rwork( * )
228  COMPLEX*16 work( * ), z( ldz, * )
229 * ..
230 *
231 * =====================================================================
232 *
233 * .. Parameters ..
234  DOUBLE PRECISION zero, one, two
235  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
236 * ..
237 * .. Local Scalars ..
238  LOGICAL lquery
239  INTEGER finish, i, icompz, ii, j, k, lgn, liwmin, ll,
240  $ lrwmin, lwmin, m, smlsiz, start
241  DOUBLE PRECISION eps, orgnrm, p, tiny
242 * ..
243 * .. External Functions ..
244  LOGICAL lsame
245  INTEGER ilaenv
246  DOUBLE PRECISION dlamch, dlanst
247  EXTERNAL lsame, ilaenv, dlamch, dlanst
248 * ..
249 * .. External Subroutines ..
250  EXTERNAL dlascl, dlaset, dstedc, dsteqr, dsterf, xerbla,
252 * ..
253 * .. Intrinsic Functions ..
254  INTRINSIC abs, dble, int, log, max, mod, sqrt
255 * ..
256 * .. Executable Statements ..
257 *
258 * Test the input parameters.
259 *
260  info = 0
261  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
262 *
263  IF( lsame( compz, 'N' ) ) THEN
264  icompz = 0
265  ELSE IF( lsame( compz, 'V' ) ) THEN
266  icompz = 1
267  ELSE IF( lsame( compz, 'I' ) ) THEN
268  icompz = 2
269  ELSE
270  icompz = -1
271  END IF
272  IF( icompz.LT.0 ) THEN
273  info = -1
274  ELSE IF( n.LT.0 ) THEN
275  info = -2
276  ELSE IF( ( ldz.LT.1 ) .OR.
277  $ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) ) THEN
278  info = -6
279  END IF
280 *
281  IF( info.EQ.0 ) THEN
282 *
283 * Compute the workspace requirements
284 *
285  smlsiz = ilaenv( 9, 'ZSTEDC', ' ', 0, 0, 0, 0 )
286  IF( n.LE.1 .OR. icompz.EQ.0 ) THEN
287  lwmin = 1
288  liwmin = 1
289  lrwmin = 1
290  ELSE IF( n.LE.smlsiz ) THEN
291  lwmin = 1
292  liwmin = 1
293  lrwmin = 2*( n - 1 )
294  ELSE IF( icompz.EQ.1 ) THEN
295  lgn = int( log( dble( n ) ) / log( two ) )
296  IF( 2**lgn.LT.n )
297  $ lgn = lgn + 1
298  IF( 2**lgn.LT.n )
299  $ lgn = lgn + 1
300  lwmin = n*n
301  lrwmin = 1 + 3*n + 2*n*lgn + 4*n**2
302  liwmin = 6 + 6*n + 5*n*lgn
303  ELSE IF( icompz.EQ.2 ) THEN
304  lwmin = 1
305  lrwmin = 1 + 4*n + 2*n**2
306  liwmin = 3 + 5*n
307  END IF
308  work( 1 ) = lwmin
309  rwork( 1 ) = lrwmin
310  iwork( 1 ) = liwmin
311 *
312  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
313  info = -8
314  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
315  info = -10
316  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
317  info = -12
318  END IF
319  END IF
320 *
321  IF( info.NE.0 ) THEN
322  CALL xerbla( 'ZSTEDC', -info )
323  return
324  ELSE IF( lquery ) THEN
325  return
326  END IF
327 *
328 * Quick return if possible
329 *
330  IF( n.EQ.0 )
331  $ return
332  IF( n.EQ.1 ) THEN
333  IF( icompz.NE.0 )
334  $ z( 1, 1 ) = one
335  return
336  END IF
337 *
338 * If the following conditional clause is removed, then the routine
339 * will use the Divide and Conquer routine to compute only the
340 * eigenvalues, which requires (3N + 3N**2) real workspace and
341 * (2 + 5N + 2N lg(N)) integer workspace.
342 * Since on many architectures DSTERF is much faster than any other
343 * algorithm for finding eigenvalues only, it is used here
344 * as the default. If the conditional clause is removed, then
345 * information on the size of workspace needs to be changed.
346 *
347 * If COMPZ = 'N', use DSTERF to compute the eigenvalues.
348 *
349  IF( icompz.EQ.0 ) THEN
350  CALL dsterf( n, d, e, info )
351  go to 70
352  END IF
353 *
354 * If N is smaller than the minimum divide size (SMLSIZ+1), then
355 * solve the problem with another solver.
356 *
357  IF( n.LE.smlsiz ) THEN
358 *
359  CALL zsteqr( compz, n, d, e, z, ldz, rwork, info )
360 *
361  ELSE
362 *
363 * If COMPZ = 'I', we simply call DSTEDC instead.
364 *
365  IF( icompz.EQ.2 ) THEN
366  CALL dlaset( 'Full', n, n, zero, one, rwork, n )
367  ll = n*n + 1
368  CALL dstedc( 'I', n, d, e, rwork, n,
369  $ rwork( ll ), lrwork-ll+1, iwork, liwork, info )
370  DO 20 j = 1, n
371  DO 10 i = 1, n
372  z( i, j ) = rwork( ( j-1 )*n+i )
373  10 continue
374  20 continue
375  go to 70
376  END IF
377 *
378 * From now on, only option left to be handled is COMPZ = 'V',
379 * i.e. ICOMPZ = 1.
380 *
381 * Scale.
382 *
383  orgnrm = dlanst( 'M', n, d, e )
384  IF( orgnrm.EQ.zero )
385  $ go to 70
386 *
387  eps = dlamch( 'Epsilon' )
388 *
389  start = 1
390 *
391 * while ( START <= N )
392 *
393  30 continue
394  IF( start.LE.n ) THEN
395 *
396 * Let FINISH be the position of the next subdiagonal entry
397 * such that E( FINISH ) <= TINY or FINISH = N if no such
398 * subdiagonal exists. The matrix identified by the elements
399 * between START and FINISH constitutes an independent
400 * sub-problem.
401 *
402  finish = start
403  40 continue
404  IF( finish.LT.n ) THEN
405  tiny = eps*sqrt( abs( d( finish ) ) )*
406  $ sqrt( abs( d( finish+1 ) ) )
407  IF( abs( e( finish ) ).GT.tiny ) THEN
408  finish = finish + 1
409  go to 40
410  END IF
411  END IF
412 *
413 * (Sub) Problem determined. Compute its size and solve it.
414 *
415  m = finish - start + 1
416  IF( m.GT.smlsiz ) THEN
417 *
418 * Scale.
419 *
420  orgnrm = dlanst( 'M', m, d( start ), e( start ) )
421  CALL dlascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
422  $ info )
423  CALL dlascl( 'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
424  $ m-1, info )
425 *
426  CALL zlaed0( n, m, d( start ), e( start ), z( 1, start ),
427  $ ldz, work, n, rwork, iwork, info )
428  IF( info.GT.0 ) THEN
429  info = ( info / ( m+1 )+start-1 )*( n+1 ) +
430  $ mod( info, ( m+1 ) ) + start - 1
431  go to 70
432  END IF
433 *
434 * Scale back.
435 *
436  CALL dlascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
437  $ info )
438 *
439  ELSE
440  CALL dsteqr( 'I', m, d( start ), e( start ), rwork, m,
441  $ rwork( m*m+1 ), info )
442  CALL zlacrm( n, m, z( 1, start ), ldz, rwork, m, work, n,
443  $ rwork( m*m+1 ) )
444  CALL zlacpy( 'A', n, m, work, n, z( 1, start ), ldz )
445  IF( info.GT.0 ) THEN
446  info = start*( n+1 ) + finish
447  go to 70
448  END IF
449  END IF
450 *
451  start = finish + 1
452  go to 30
453  END IF
454 *
455 * endwhile
456 *
457 * If the problem split any number of times, then the eigenvalues
458 * will not be properly ordered. Here we permute the eigenvalues
459 * (and the associated eigenvectors) into ascending order.
460 *
461  IF( m.NE.n ) THEN
462 *
463 * Use Selection Sort to minimize swaps of eigenvectors
464 *
465  DO 60 ii = 2, n
466  i = ii - 1
467  k = i
468  p = d( i )
469  DO 50 j = ii, n
470  IF( d( j ).LT.p ) THEN
471  k = j
472  p = d( j )
473  END IF
474  50 continue
475  IF( k.NE.i ) THEN
476  d( k ) = d( i )
477  d( i ) = p
478  CALL zswap( n, z( 1, i ), 1, z( 1, k ), 1 )
479  END IF
480  60 continue
481  END IF
482  END IF
483 *
484  70 continue
485  work( 1 ) = lwmin
486  rwork( 1 ) = lrwmin
487  iwork( 1 ) = liwmin
488 *
489  return
490 *
491 * End of ZSTEDC
492 *
493  END