LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
claed0.f
Go to the documentation of this file.
1 *> \brief \b CLAED0 used by CSTEDC. 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 CLAED0 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claed0.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claed0.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claed0.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CLAED0( 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 * REAL D( * ), E( * ), RWORK( * )
30 * COMPLEX Q( LDQ, * ), QSTORE( LDQS, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> Using the divide and conquer method, CLAED0 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 REAL 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 REAL 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 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 REAL 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 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 *> \ingroup complexOTHERcomputational
141 *
142 * =====================================================================
143  SUBROUTINE claed0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
144  $ IWORK, INFO )
145 *
146 * -- LAPACK computational routine --
147 * -- LAPACK is a software package provided by Univ. of Tennessee, --
148 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
149 *
150 * .. Scalar Arguments ..
151  INTEGER INFO, LDQ, LDQS, N, QSIZ
152 * ..
153 * .. Array Arguments ..
154  INTEGER IWORK( * )
155  REAL D( * ), E( * ), RWORK( * )
156  COMPLEX Q( LDQ, * ), QSTORE( LDQS, * )
157 * ..
158 *
159 * =====================================================================
160 *
161 * Warning: N could be as big as QSIZ!
162 *
163 * .. Parameters ..
164  REAL TWO
165  parameter( two = 2.e+0 )
166 * ..
167 * .. Local Scalars ..
168  INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
169  $ igivpt, indxq, iperm, iprmpt, iq, iqptr, iwrem,
170  $ j, k, lgn, ll, matsiz, msd2, smlsiz, smm1,
171  $ spm1, spm2, submat, subpbs, tlvls
172  REAL TEMP
173 * ..
174 * .. External Subroutines ..
175  EXTERNAL ccopy, clacrm, claed7, scopy, ssteqr, xerbla
176 * ..
177 * .. External Functions ..
178  INTEGER ILAENV
179  EXTERNAL ilaenv
180 * ..
181 * .. Intrinsic Functions ..
182  INTRINSIC abs, int, log, max, real
183 * ..
184 * .. Executable Statements ..
185 *
186 * Test the input parameters.
187 *
188  info = 0
189 *
190 * IF( ICOMPQ .LT. 0 .OR. ICOMPQ .GT. 2 ) THEN
191 * INFO = -1
192 * ELSE IF( ( ICOMPQ .EQ. 1 ) .AND. ( QSIZ .LT. MAX( 0, N ) ) )
193 * $ THEN
194  IF( qsiz.LT.max( 0, n ) ) THEN
195  info = -1
196  ELSE IF( n.LT.0 ) THEN
197  info = -2
198  ELSE IF( ldq.LT.max( 1, n ) ) THEN
199  info = -6
200  ELSE IF( ldqs.LT.max( 1, n ) ) THEN
201  info = -8
202  END IF
203  IF( info.NE.0 ) THEN
204  CALL xerbla( 'CLAED0', -info )
205  RETURN
206  END IF
207 *
208 * Quick return if possible
209 *
210  IF( n.EQ.0 )
211  $ RETURN
212 *
213  smlsiz = ilaenv( 9, 'CLAED0', ' ', 0, 0, 0, 0 )
214 *
215 * Determine the size and placement of the submatrices, and save in
216 * the leading elements of IWORK.
217 *
218  iwork( 1 ) = n
219  subpbs = 1
220  tlvls = 0
221  10 CONTINUE
222  IF( iwork( subpbs ).GT.smlsiz ) THEN
223  DO 20 j = subpbs, 1, -1
224  iwork( 2*j ) = ( iwork( j )+1 ) / 2
225  iwork( 2*j-1 ) = iwork( j ) / 2
226  20 CONTINUE
227  tlvls = tlvls + 1
228  subpbs = 2*subpbs
229  GO TO 10
230  END IF
231  DO 30 j = 2, subpbs
232  iwork( j ) = iwork( j ) + iwork( j-1 )
233  30 CONTINUE
234 *
235 * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
236 * using rank-1 modifications (cuts).
237 *
238  spm1 = subpbs - 1
239  DO 40 i = 1, spm1
240  submat = iwork( i ) + 1
241  smm1 = submat - 1
242  d( smm1 ) = d( smm1 ) - abs( e( smm1 ) )
243  d( submat ) = d( submat ) - abs( e( smm1 ) )
244  40 CONTINUE
245 *
246  indxq = 4*n + 3
247 *
248 * Set up workspaces for eigenvalues only/accumulate new vectors
249 * routine
250 *
251  temp = log( real( n ) ) / log( two )
252  lgn = int( temp )
253  IF( 2**lgn.LT.n )
254  $ lgn = lgn + 1
255  IF( 2**lgn.LT.n )
256  $ lgn = lgn + 1
257  iprmpt = indxq + n + 1
258  iperm = iprmpt + n*lgn
259  iqptr = iperm + n*lgn
260  igivpt = iqptr + n + 2
261  igivcl = igivpt + n*lgn
262 *
263  igivnm = 1
264  iq = igivnm + 2*n*lgn
265  iwrem = iq + n**2 + 1
266 * Initialize pointers
267  DO 50 i = 0, subpbs
268  iwork( iprmpt+i ) = 1
269  iwork( igivpt+i ) = 1
270  50 CONTINUE
271  iwork( iqptr ) = 1
272 *
273 * Solve each submatrix eigenproblem at the bottom of the divide and
274 * conquer tree.
275 *
276  curr = 0
277  DO 70 i = 0, spm1
278  IF( i.EQ.0 ) THEN
279  submat = 1
280  matsiz = iwork( 1 )
281  ELSE
282  submat = iwork( i ) + 1
283  matsiz = iwork( i+1 ) - iwork( i )
284  END IF
285  ll = iq - 1 + iwork( iqptr+curr )
286  CALL ssteqr( 'I', matsiz, d( submat ), e( submat ),
287  $ rwork( ll ), matsiz, rwork, info )
288  CALL clacrm( qsiz, matsiz, q( 1, submat ), ldq, rwork( ll ),
289  $ matsiz, qstore( 1, submat ), ldqs,
290  $ rwork( iwrem ) )
291  iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2
292  curr = curr + 1
293  IF( info.GT.0 ) THEN
294  info = submat*( n+1 ) + submat + matsiz - 1
295  RETURN
296  END IF
297  k = 1
298  DO 60 j = submat, iwork( i+1 )
299  iwork( indxq+j ) = k
300  k = k + 1
301  60 CONTINUE
302  70 CONTINUE
303 *
304 * Successively merge eigensystems of adjacent submatrices
305 * into eigensystem for the corresponding larger matrix.
306 *
307 * while ( SUBPBS > 1 )
308 *
309  curlvl = 1
310  80 CONTINUE
311  IF( subpbs.GT.1 ) THEN
312  spm2 = subpbs - 2
313  DO 90 i = 0, spm2, 2
314  IF( i.EQ.0 ) THEN
315  submat = 1
316  matsiz = iwork( 2 )
317  msd2 = iwork( 1 )
318  curprb = 0
319  ELSE
320  submat = iwork( i ) + 1
321  matsiz = iwork( i+2 ) - iwork( i )
322  msd2 = matsiz / 2
323  curprb = curprb + 1
324  END IF
325 *
326 * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
327 * into an eigensystem of size MATSIZ. CLAED7 handles the case
328 * when the eigenvectors of a full or band Hermitian matrix (which
329 * was reduced to tridiagonal form) are desired.
330 *
331 * I am free to use Q as a valuable working space until Loop 150.
332 *
333  CALL claed7( matsiz, msd2, qsiz, tlvls, curlvl, curprb,
334  $ d( submat ), qstore( 1, submat ), ldqs,
335  $ e( submat+msd2-1 ), iwork( indxq+submat ),
336  $ rwork( iq ), iwork( iqptr ), iwork( iprmpt ),
337  $ iwork( iperm ), iwork( igivpt ),
338  $ iwork( igivcl ), rwork( igivnm ),
339  $ q( 1, submat ), rwork( iwrem ),
340  $ iwork( subpbs+1 ), info )
341  IF( info.GT.0 ) THEN
342  info = submat*( n+1 ) + submat + matsiz - 1
343  RETURN
344  END IF
345  iwork( i / 2+1 ) = iwork( i+2 )
346  90 CONTINUE
347  subpbs = subpbs / 2
348  curlvl = curlvl + 1
349  GO TO 80
350  END IF
351 *
352 * end while
353 *
354 * Re-merge the eigenvalues/vectors which were deflated at the final
355 * merge step.
356 *
357  DO 100 i = 1, n
358  j = iwork( indxq+i )
359  rwork( i ) = d( j )
360  CALL ccopy( qsiz, qstore( 1, j ), 1, q( 1, i ), 1 )
361  100 CONTINUE
362  CALL scopy( n, rwork, 1, d, 1 )
363 *
364  RETURN
365 *
366 * End of CLAED0
367 *
368  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:131
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine clacrm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
CLACRM multiplies a complex matrix by a square real matrix.
Definition: clacrm.f:114
subroutine claed0(QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK, IWORK, INFO)
CLAED0 used by CSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
Definition: claed0.f:145
subroutine claed7(N, CUTPNT, QSIZ, TLVLS, CURLVL, CURPBM, D, Q, LDQ, RHO, INDXQ, QSTORE, QPTR, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, WORK, RWORK, IWORK, INFO)
CLAED7 used by CSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a ...
Definition: claed7.f:249
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82