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