LAPACK  3.8.0 LAPACK: Linear Algebra PACKage
claed0.f
Go to the documentation of this file.
1 *> \brief \b CLAED0 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
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 *> \date December 2016
141 *
142 *> \ingroup complexOTHERcomputational
143 *
144 * =====================================================================
145  SUBROUTINE claed0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
146  \$ IWORK, INFO )
147 *
148 * -- LAPACK computational routine (version 3.7.0) --
149 * -- LAPACK is a software package provided by Univ. of Tennessee, --
150 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
151 * December 2016
152 *
153 * .. Scalar Arguments ..
154  INTEGER INFO, LDQ, LDQS, N, QSIZ
155 * ..
156 * .. Array Arguments ..
157  INTEGER IWORK( * )
158  REAL D( * ), E( * ), RWORK( * )
159  COMPLEX Q( ldq, * ), QSTORE( ldqs, * )
160 * ..
161 *
162 * =====================================================================
163 *
164 * Warning: N could be as big as QSIZ!
165 *
166 * .. Parameters ..
167  REAL TWO
168  parameter( two = 2.e+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  REAL TEMP
176 * ..
177 * .. External Subroutines ..
178  EXTERNAL ccopy, clacrm, claed7, scopy, ssteqr, xerbla
179 * ..
180 * .. External Functions ..
181  INTEGER ILAENV
182  EXTERNAL ilaenv
183 * ..
184 * .. Intrinsic Functions ..
185  INTRINSIC abs, int, log, max, real
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( 'CLAED0', -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, 'CLAED0', ' ', 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( REAL( 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 ssteqr( 'I', matsiz, d( submat ), e( submat ),
290  \$ rwork( ll ), matsiz, rwork, info )
291  CALL clacrm( 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. CLAED7 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 claed7( 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 ccopy( qsiz, qstore( 1, j ), 1, q( 1, i ), 1 )
364  100 CONTINUE
365  CALL scopy( n, rwork, 1, d, 1 )
366 *
367  RETURN
368 *
369 * End of CLAED0
370 *
371  END
subroutine claed0(QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK, IWORK, INFO)
CLAED0 used by sstedc. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
Definition: claed0.f:147
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84
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 sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a ...
Definition: claed7.f:251
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:133
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clacrm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
CLACRM multiplies a complex matrix by a square real matrix.
Definition: clacrm.f:116