LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dstedc.f
Go to the documentation of this file.
1 *> \brief \b DSTEDC
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, 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 N <= 1 then LWORK must be at least 1.
117 *> If COMPZ = 'V' and N > 1 then LWORK must be at least
118 *> ( 1 + 3*N + 2*N*lg N + 4*N**2 ),
119 *> where lg( N ) = smallest integer k such
120 *> that 2**k >= N.
121 *> If COMPZ = 'I' and N > 1 then LWORK must be at least
122 *> ( 1 + 4*N + N**2 ).
123 *> Note that for COMPZ = 'I' or 'V', then if N is less than or
124 *> equal to the minimum divide size, usually 25, then LWORK need
125 *> only be max(1,2*(N-1)).
126 *>
127 *> If LWORK = -1, then a workspace query is assumed; the routine
128 *> only calculates the optimal size of the WORK array, returns
129 *> this value as the first entry of the WORK array, and no error
130 *> message related to LWORK is issued by XERBLA.
131 *> \endverbatim
132 *>
133 *> \param[out] IWORK
134 *> \verbatim
135 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
136 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
137 *> \endverbatim
138 *>
139 *> \param[in] LIWORK
140 *> \verbatim
141 *> LIWORK is INTEGER
142 *> The dimension of the array IWORK.
143 *> If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
144 *> If COMPZ = 'V' and N > 1 then LIWORK must be at least
145 *> ( 6 + 6*N + 5*N*lg N ).
146 *> If COMPZ = 'I' and N > 1 then LIWORK must be at least
147 *> ( 3 + 5*N ).
148 *> Note that for COMPZ = 'I' or 'V', then if N is less than or
149 *> equal to the minimum divide size, usually 25, then LIWORK
150 *> need only be 1.
151 *>
152 *> If LIWORK = -1, then a workspace query is assumed; the
153 *> routine only calculates the optimal size of the IWORK array,
154 *> returns this value as the first entry of the IWORK array, and
155 *> no error message related to LIWORK is issued by XERBLA.
156 *> \endverbatim
157 *>
158 *> \param[out] INFO
159 *> \verbatim
160 *> INFO is INTEGER
161 *> = 0: successful exit.
162 *> < 0: if INFO = -i, the i-th argument had an illegal value.
163 *> > 0: The algorithm failed to compute an eigenvalue while
164 *> working on the submatrix lying in rows and columns
165 *> INFO/(N+1) through mod(INFO,N+1).
166 *> \endverbatim
167 *
168 * Authors:
169 * ========
170 *
171 *> \author Univ. of Tennessee
172 *> \author Univ. of California Berkeley
173 *> \author Univ. of Colorado Denver
174 *> \author NAG Ltd.
175 *
176 *> \ingroup auxOTHERcomputational
177 *
178 *> \par Contributors:
179 * ==================
180 *>
181 *> Jeff Rutter, Computer Science Division, University of California
182 *> at Berkeley, USA \n
183 *> Modified by Francoise Tisseur, University of Tennessee
184 *>
185 * =====================================================================
186  SUBROUTINE dstedc( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
187  $ LIWORK, INFO )
188 *
189 * -- LAPACK computational routine --
190 * -- LAPACK is a software package provided by Univ. of Tennessee, --
191 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
192 *
193 * .. Scalar Arguments ..
194  CHARACTER COMPZ
195  INTEGER INFO, LDZ, LIWORK, LWORK, N
196 * ..
197 * .. Array Arguments ..
198  INTEGER IWORK( * )
199  DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
200 * ..
201 *
202 * =====================================================================
203 *
204 * .. Parameters ..
205  DOUBLE PRECISION ZERO, ONE, TWO
206  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
207 * ..
208 * .. Local Scalars ..
209  LOGICAL LQUERY
210  INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
211  $ lwmin, m, smlsiz, start, storez, strtrw
212  DOUBLE PRECISION EPS, ORGNRM, P, TINY
213 * ..
214 * .. External Functions ..
215  LOGICAL LSAME
216  INTEGER ILAENV
217  DOUBLE PRECISION DLAMCH, DLANST
218  EXTERNAL lsame, ilaenv, dlamch, dlanst
219 * ..
220 * .. External Subroutines ..
221  EXTERNAL dgemm, dlacpy, dlaed0, dlascl, dlaset, dlasrt,
223 * ..
224 * .. Intrinsic Functions ..
225  INTRINSIC abs, dble, int, log, max, mod, sqrt
226 * ..
227 * .. Executable Statements ..
228 *
229 * Test the input parameters.
230 *
231  info = 0
232  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
233 *
234  IF( lsame( compz, 'N' ) ) THEN
235  icompz = 0
236  ELSE IF( lsame( compz, 'V' ) ) THEN
237  icompz = 1
238  ELSE IF( lsame( compz, 'I' ) ) THEN
239  icompz = 2
240  ELSE
241  icompz = -1
242  END IF
243  IF( icompz.LT.0 ) THEN
244  info = -1
245  ELSE IF( n.LT.0 ) THEN
246  info = -2
247  ELSE IF( ( ldz.LT.1 ) .OR.
248  $ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) ) THEN
249  info = -6
250  END IF
251 *
252  IF( info.EQ.0 ) THEN
253 *
254 * Compute the workspace requirements
255 *
256  smlsiz = ilaenv( 9, 'DSTEDC', ' ', 0, 0, 0, 0 )
257  IF( n.LE.1 .OR. icompz.EQ.0 ) THEN
258  liwmin = 1
259  lwmin = 1
260  ELSE IF( n.LE.smlsiz ) THEN
261  liwmin = 1
262  lwmin = 2*( n - 1 )
263  ELSE
264  lgn = int( log( dble( n ) )/log( two ) )
265  IF( 2**lgn.LT.n )
266  $ lgn = lgn + 1
267  IF( 2**lgn.LT.n )
268  $ lgn = lgn + 1
269  IF( icompz.EQ.1 ) THEN
270  lwmin = 1 + 3*n + 2*n*lgn + 4*n**2
271  liwmin = 6 + 6*n + 5*n*lgn
272  ELSE IF( icompz.EQ.2 ) THEN
273  lwmin = 1 + 4*n + n**2
274  liwmin = 3 + 5*n
275  END IF
276  END IF
277  work( 1 ) = lwmin
278  iwork( 1 ) = liwmin
279 *
280  IF( lwork.LT.lwmin .AND. .NOT. lquery ) THEN
281  info = -8
282  ELSE IF( liwork.LT.liwmin .AND. .NOT. lquery ) THEN
283  info = -10
284  END IF
285  END IF
286 *
287  IF( info.NE.0 ) THEN
288  CALL xerbla( 'DSTEDC', -info )
289  RETURN
290  ELSE IF (lquery) THEN
291  RETURN
292  END IF
293 *
294 * Quick return if possible
295 *
296  IF( n.EQ.0 )
297  $ RETURN
298  IF( n.EQ.1 ) THEN
299  IF( icompz.NE.0 )
300  $ z( 1, 1 ) = one
301  RETURN
302  END IF
303 *
304 * If the following conditional clause is removed, then the routine
305 * will use the Divide and Conquer routine to compute only the
306 * eigenvalues, which requires (3N + 3N**2) real workspace and
307 * (2 + 5N + 2N lg(N)) integer workspace.
308 * Since on many architectures DSTERF is much faster than any other
309 * algorithm for finding eigenvalues only, it is used here
310 * as the default. If the conditional clause is removed, then
311 * information on the size of workspace needs to be changed.
312 *
313 * If COMPZ = 'N', use DSTERF to compute the eigenvalues.
314 *
315  IF( icompz.EQ.0 ) THEN
316  CALL dsterf( n, d, e, info )
317  GO TO 50
318  END IF
319 *
320 * If N is smaller than the minimum divide size (SMLSIZ+1), then
321 * solve the problem with another solver.
322 *
323  IF( n.LE.smlsiz ) THEN
324 *
325  CALL dsteqr( compz, n, d, e, z, ldz, work, info )
326 *
327  ELSE
328 *
329 * If COMPZ = 'V', the Z matrix must be stored elsewhere for later
330 * use.
331 *
332  IF( icompz.EQ.1 ) THEN
333  storez = 1 + n*n
334  ELSE
335  storez = 1
336  END IF
337 *
338  IF( icompz.EQ.2 ) THEN
339  CALL dlaset( 'Full', n, n, zero, one, z, ldz )
340  END IF
341 *
342 * Scale.
343 *
344  orgnrm = dlanst( 'M', n, d, e )
345  IF( orgnrm.EQ.zero )
346  $ GO TO 50
347 *
348  eps = dlamch( 'Epsilon' )
349 *
350  start = 1
351 *
352 * while ( START <= N )
353 *
354  10 CONTINUE
355  IF( start.LE.n ) THEN
356 *
357 * Let FINISH be the position of the next subdiagonal entry
358 * such that E( FINISH ) <= TINY or FINISH = N if no such
359 * subdiagonal exists. The matrix identified by the elements
360 * between START and FINISH constitutes an independent
361 * sub-problem.
362 *
363  finish = start
364  20 CONTINUE
365  IF( finish.LT.n ) THEN
366  tiny = eps*sqrt( abs( d( finish ) ) )*
367  $ sqrt( abs( d( finish+1 ) ) )
368  IF( abs( e( finish ) ).GT.tiny ) THEN
369  finish = finish + 1
370  GO TO 20
371  END IF
372  END IF
373 *
374 * (Sub) Problem determined. Compute its size and solve it.
375 *
376  m = finish - start + 1
377  IF( m.EQ.1 ) THEN
378  start = finish + 1
379  GO TO 10
380  END IF
381  IF( m.GT.smlsiz ) THEN
382 *
383 * Scale.
384 *
385  orgnrm = dlanst( 'M', m, d( start ), e( start ) )
386  CALL dlascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
387  $ info )
388  CALL dlascl( 'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
389  $ m-1, info )
390 *
391  IF( icompz.EQ.1 ) THEN
392  strtrw = 1
393  ELSE
394  strtrw = start
395  END IF
396  CALL dlaed0( icompz, n, m, d( start ), e( start ),
397  $ z( strtrw, start ), ldz, work( 1 ), n,
398  $ work( storez ), iwork, info )
399  IF( info.NE.0 ) THEN
400  info = ( info / ( m+1 )+start-1 )*( n+1 ) +
401  $ mod( info, ( m+1 ) ) + start - 1
402  GO TO 50
403  END IF
404 *
405 * Scale back.
406 *
407  CALL dlascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
408  $ info )
409 *
410  ELSE
411  IF( icompz.EQ.1 ) THEN
412 *
413 * Since QR won't update a Z matrix which is larger than
414 * the length of D, we must solve the sub-problem in a
415 * workspace and then multiply back into Z.
416 *
417  CALL dsteqr( 'I', m, d( start ), e( start ), work, m,
418  $ work( m*m+1 ), info )
419  CALL dlacpy( 'A', n, m, z( 1, start ), ldz,
420  $ work( storez ), n )
421  CALL dgemm( 'N', 'N', n, m, m, one,
422  $ work( storez ), n, work, m, zero,
423  $ z( 1, start ), ldz )
424  ELSE IF( icompz.EQ.2 ) THEN
425  CALL dsteqr( 'I', m, d( start ), e( start ),
426  $ z( start, start ), ldz, work, info )
427  ELSE
428  CALL dsterf( m, d( start ), e( start ), info )
429  END IF
430  IF( info.NE.0 ) THEN
431  info = start*( n+1 ) + finish
432  GO TO 50
433  END IF
434  END IF
435 *
436  start = finish + 1
437  GO TO 10
438  END IF
439 *
440 * endwhile
441 *
442  IF( icompz.EQ.0 ) THEN
443 *
444 * Use Quick Sort
445 *
446  CALL dlasrt( 'I', n, d, info )
447 *
448  ELSE
449 *
450 * Use Selection Sort to minimize swaps of eigenvectors
451 *
452  DO 40 ii = 2, n
453  i = ii - 1
454  k = i
455  p = d( i )
456  DO 30 j = ii, n
457  IF( d( j ).LT.p ) THEN
458  k = j
459  p = d( j )
460  END IF
461  30 CONTINUE
462  IF( k.NE.i ) THEN
463  d( k ) = d( i )
464  d( i ) = p
465  CALL dswap( n, z( 1, i ), 1, z( 1, k ), 1 )
466  END IF
467  40 CONTINUE
468  END IF
469  END IF
470 *
471  50 CONTINUE
472  work( 1 ) = lwmin
473  iwork( 1 ) = liwmin
474 *
475  RETURN
476 *
477 * End of DSTEDC
478 *
479  END
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:143
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:131
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:88
subroutine dstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
DSTEDC
Definition: dstedc.f:188
subroutine dlaed0(ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, WORK, IWORK, INFO)
DLAED0 used by DSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
Definition: dlaed0.f:172
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187