LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zlaed0.f
Go to the documentation of this file.
1 *> \brief \b ZLAED0 used by sstedc. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZLAED0 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaed0.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaed0.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaed0.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
22 * IWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER INFO, LDQ, LDQS, N, QSIZ
26 * ..
27 * .. Array Arguments ..
28 * INTEGER IWORK( * )
29 * DOUBLE PRECISION D( * ), E( * ), RWORK( * )
30 * COMPLEX*16 Q( LDQ, * ), QSTORE( LDQS, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> Using the divide and conquer method, ZLAED0 computes all eigenvalues
40 *> of a symmetric tridiagonal matrix which is one diagonal block of
41 *> those from reducing a dense or band Hermitian matrix and
42 *> corresponding eigenvectors of the dense or band matrix.
43 *> \endverbatim
44 *
45 * Arguments:
46 * ==========
47 *
48 *> \param[in] QSIZ
49 *> \verbatim
50 *> QSIZ is INTEGER
51 *> The dimension of the unitary matrix used to reduce
52 *> the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
53 *> \endverbatim
54 *>
55 *> \param[in] N
56 *> \verbatim
57 *> N is INTEGER
58 *> The dimension of the symmetric tridiagonal matrix. N >= 0.
59 *> \endverbatim
60 *>
61 *> \param[in,out] D
62 *> \verbatim
63 *> D is DOUBLE PRECISION array, dimension (N)
64 *> On entry, the diagonal elements of the tridiagonal matrix.
65 *> On exit, the eigenvalues in ascending order.
66 *> \endverbatim
67 *>
68 *> \param[in,out] E
69 *> \verbatim
70 *> E is DOUBLE PRECISION array, dimension (N-1)
71 *> On entry, the off-diagonal elements of the tridiagonal matrix.
72 *> On exit, E has been destroyed.
73 *> \endverbatim
74 *>
75 *> \param[in,out] Q
76 *> \verbatim
77 *> Q is COMPLEX*16 array, dimension (LDQ,N)
78 *> On entry, Q must contain an QSIZ x N matrix whose columns
79 *> unitarily orthonormal. It is a part of the unitary matrix
80 *> that reduces the full dense Hermitian matrix to a
81 *> (reducible) symmetric tridiagonal matrix.
82 *> \endverbatim
83 *>
84 *> \param[in] LDQ
85 *> \verbatim
86 *> LDQ is INTEGER
87 *> The leading dimension of the array Q. LDQ >= max(1,N).
88 *> \endverbatim
89 *>
90 *> \param[out] IWORK
91 *> \verbatim
92 *> IWORK is INTEGER array,
93 *> the dimension of IWORK must be at least
94 *> 6 + 6*N + 5*N*lg N
95 *> ( lg( N ) = smallest integer k
96 *> such that 2^k >= N )
97 *> \endverbatim
98 *>
99 *> \param[out] RWORK
100 *> \verbatim
101 *> RWORK is DOUBLE PRECISION array,
102 *> dimension (1 + 3*N + 2*N*lg N + 3*N**2)
103 *> ( lg( N ) = smallest integer k
104 *> such that 2^k >= N )
105 *> \endverbatim
106 *>
107 *> \param[out] QSTORE
108 *> \verbatim
109 *> QSTORE is COMPLEX*16 array, dimension (LDQS, N)
110 *> Used to store parts of
111 *> the eigenvector matrix when the updating matrix multiplies
112 *> take place.
113 *> \endverbatim
114 *>
115 *> \param[in] LDQS
116 *> \verbatim
117 *> LDQS is INTEGER
118 *> The leading dimension of the array QSTORE.
119 *> LDQS >= max(1,N).
120 *> \endverbatim
121 *>
122 *> \param[out] INFO
123 *> \verbatim
124 *> INFO is INTEGER
125 *> = 0: successful exit.
126 *> < 0: if INFO = -i, the i-th argument had an illegal value.
127 *> > 0: The algorithm failed to compute an eigenvalue while
128 *> working on the submatrix lying in rows and columns
129 *> INFO/(N+1) through mod(INFO,N+1).
130 *> \endverbatim
131 *
132 * Authors:
133 * ========
134 *
135 *> \author Univ. of Tennessee
136 *> \author Univ. of California Berkeley
137 *> \author Univ. of Colorado Denver
138 *> \author NAG Ltd.
139 *
140 *> \date September 2012
141 *
142 *> \ingroup complex16OTHERcomputational
143 *
144 * =====================================================================
145  SUBROUTINE zlaed0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
146  $ iwork, info )
147 *
148 * -- LAPACK computational routine (version 3.4.2) --
149 * -- LAPACK is a software package provided by Univ. of Tennessee, --
150 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
151 * September 2012
152 *
153 * .. Scalar Arguments ..
154  INTEGER info, ldq, ldqs, n, qsiz
155 * ..
156 * .. Array Arguments ..
157  INTEGER iwork( * )
158  DOUBLE PRECISION d( * ), e( * ), rwork( * )
159  COMPLEX*16 q( ldq, * ), qstore( ldqs, * )
160 * ..
161 *
162 * =====================================================================
163 *
164 * Warning: N could be as big as QSIZ!
165 *
166 * .. Parameters ..
167  DOUBLE PRECISION two
168  parameter( two = 2.d+0 )
169 * ..
170 * .. Local Scalars ..
171  INTEGER curlvl, curprb, curr, i, igivcl, igivnm,
172  $ igivpt, indxq, iperm, iprmpt, iq, iqptr, iwrem,
173  $ j, k, lgn, ll, matsiz, msd2, smlsiz, smm1,
174  $ spm1, spm2, submat, subpbs, tlvls
175  DOUBLE PRECISION temp
176 * ..
177 * .. External Subroutines ..
178  EXTERNAL dcopy, dsteqr, xerbla, zcopy, zlacrm, zlaed7
179 * ..
180 * .. External Functions ..
181  INTEGER ilaenv
182  EXTERNAL ilaenv
183 * ..
184 * .. Intrinsic Functions ..
185  INTRINSIC abs, dble, int, log, max
186 * ..
187 * .. Executable Statements ..
188 *
189 * Test the input parameters.
190 *
191  info = 0
192 *
193 * IF( ICOMPQ .LT. 0 .OR. ICOMPQ .GT. 2 ) THEN
194 * INFO = -1
195 * ELSE IF( ( ICOMPQ .EQ. 1 ) .AND. ( QSIZ .LT. MAX( 0, N ) ) )
196 * $ THEN
197  IF( qsiz.LT.max( 0, n ) ) THEN
198  info = -1
199  ELSE IF( n.LT.0 ) THEN
200  info = -2
201  ELSE IF( ldq.LT.max( 1, n ) ) THEN
202  info = -6
203  ELSE IF( ldqs.LT.max( 1, n ) ) THEN
204  info = -8
205  END IF
206  IF( info.NE.0 ) THEN
207  CALL xerbla( 'ZLAED0', -info )
208  return
209  END IF
210 *
211 * Quick return if possible
212 *
213  IF( n.EQ.0 )
214  $ return
215 *
216  smlsiz = ilaenv( 9, 'ZLAED0', ' ', 0, 0, 0, 0 )
217 *
218 * Determine the size and placement of the submatrices, and save in
219 * the leading elements of IWORK.
220 *
221  iwork( 1 ) = n
222  subpbs = 1
223  tlvls = 0
224  10 continue
225  IF( iwork( subpbs ).GT.smlsiz ) THEN
226  DO 20 j = subpbs, 1, -1
227  iwork( 2*j ) = ( iwork( j )+1 ) / 2
228  iwork( 2*j-1 ) = iwork( j ) / 2
229  20 continue
230  tlvls = tlvls + 1
231  subpbs = 2*subpbs
232  go to 10
233  END IF
234  DO 30 j = 2, subpbs
235  iwork( j ) = iwork( j ) + iwork( j-1 )
236  30 continue
237 *
238 * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
239 * using rank-1 modifications (cuts).
240 *
241  spm1 = subpbs - 1
242  DO 40 i = 1, spm1
243  submat = iwork( i ) + 1
244  smm1 = submat - 1
245  d( smm1 ) = d( smm1 ) - abs( e( smm1 ) )
246  d( submat ) = d( submat ) - abs( e( smm1 ) )
247  40 continue
248 *
249  indxq = 4*n + 3
250 *
251 * Set up workspaces for eigenvalues only/accumulate new vectors
252 * routine
253 *
254  temp = log( dble( n ) ) / log( two )
255  lgn = int( temp )
256  IF( 2**lgn.LT.n )
257  $ lgn = lgn + 1
258  IF( 2**lgn.LT.n )
259  $ lgn = lgn + 1
260  iprmpt = indxq + n + 1
261  iperm = iprmpt + n*lgn
262  iqptr = iperm + n*lgn
263  igivpt = iqptr + n + 2
264  igivcl = igivpt + n*lgn
265 *
266  igivnm = 1
267  iq = igivnm + 2*n*lgn
268  iwrem = iq + n**2 + 1
269 * Initialize pointers
270  DO 50 i = 0, subpbs
271  iwork( iprmpt+i ) = 1
272  iwork( igivpt+i ) = 1
273  50 continue
274  iwork( iqptr ) = 1
275 *
276 * Solve each submatrix eigenproblem at the bottom of the divide and
277 * conquer tree.
278 *
279  curr = 0
280  DO 70 i = 0, spm1
281  IF( i.EQ.0 ) THEN
282  submat = 1
283  matsiz = iwork( 1 )
284  ELSE
285  submat = iwork( i ) + 1
286  matsiz = iwork( i+1 ) - iwork( i )
287  END IF
288  ll = iq - 1 + iwork( iqptr+curr )
289  CALL dsteqr( 'I', matsiz, d( submat ), e( submat ),
290  $ rwork( ll ), matsiz, rwork, info )
291  CALL zlacrm( qsiz, matsiz, q( 1, submat ), ldq, rwork( ll ),
292  $ matsiz, qstore( 1, submat ), ldqs,
293  $ rwork( iwrem ) )
294  iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2
295  curr = curr + 1
296  IF( info.GT.0 ) THEN
297  info = submat*( n+1 ) + submat + matsiz - 1
298  return
299  END IF
300  k = 1
301  DO 60 j = submat, iwork( i+1 )
302  iwork( indxq+j ) = k
303  k = k + 1
304  60 continue
305  70 continue
306 *
307 * Successively merge eigensystems of adjacent submatrices
308 * into eigensystem for the corresponding larger matrix.
309 *
310 * while ( SUBPBS > 1 )
311 *
312  curlvl = 1
313  80 continue
314  IF( subpbs.GT.1 ) THEN
315  spm2 = subpbs - 2
316  DO 90 i = 0, spm2, 2
317  IF( i.EQ.0 ) THEN
318  submat = 1
319  matsiz = iwork( 2 )
320  msd2 = iwork( 1 )
321  curprb = 0
322  ELSE
323  submat = iwork( i ) + 1
324  matsiz = iwork( i+2 ) - iwork( i )
325  msd2 = matsiz / 2
326  curprb = curprb + 1
327  END IF
328 *
329 * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
330 * into an eigensystem of size MATSIZ. ZLAED7 handles the case
331 * when the eigenvectors of a full or band Hermitian matrix (which
332 * was reduced to tridiagonal form) are desired.
333 *
334 * I am free to use Q as a valuable working space until Loop 150.
335 *
336  CALL zlaed7( matsiz, msd2, qsiz, tlvls, curlvl, curprb,
337  $ d( submat ), qstore( 1, submat ), ldqs,
338  $ e( submat+msd2-1 ), iwork( indxq+submat ),
339  $ rwork( iq ), iwork( iqptr ), iwork( iprmpt ),
340  $ iwork( iperm ), iwork( igivpt ),
341  $ iwork( igivcl ), rwork( igivnm ),
342  $ q( 1, submat ), rwork( iwrem ),
343  $ iwork( subpbs+1 ), info )
344  IF( info.GT.0 ) THEN
345  info = submat*( n+1 ) + submat + matsiz - 1
346  return
347  END IF
348  iwork( i / 2+1 ) = iwork( i+2 )
349  90 continue
350  subpbs = subpbs / 2
351  curlvl = curlvl + 1
352  go to 80
353  END IF
354 *
355 * end while
356 *
357 * Re-merge the eigenvalues/vectors which were deflated at the final
358 * merge step.
359 *
360  DO 100 i = 1, n
361  j = iwork( indxq+i )
362  rwork( i ) = d( j )
363  CALL zcopy( qsiz, qstore( 1, j ), 1, q( 1, i ), 1 )
364  100 continue
365  CALL dcopy( n, rwork, 1, d, 1 )
366 *
367  return
368 *
369 * End of ZLAED0
370 *
371  END