LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dlaed0.f
Go to the documentation of this file.
1 *> \brief \b DLAED0 used by DSTEDC. 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 DLAED0 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed0.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed0.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed0.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
22 * WORK, IWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
26 * ..
27 * .. Array Arguments ..
28 * INTEGER IWORK( * )
29 * DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
30 * $ WORK( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> DLAED0 computes all eigenvalues and corresponding eigenvectors of a
40 *> symmetric tridiagonal matrix using the divide and conquer method.
41 *> \endverbatim
42 *
43 * Arguments:
44 * ==========
45 *
46 *> \param[in] ICOMPQ
47 *> \verbatim
48 *> ICOMPQ is INTEGER
49 *> = 0: Compute eigenvalues only.
50 *> = 1: Compute eigenvectors of original dense symmetric matrix
51 *> also. On entry, Q contains the orthogonal matrix used
52 *> to reduce the original matrix to tridiagonal form.
53 *> = 2: Compute eigenvalues and eigenvectors of tridiagonal
54 *> matrix.
55 *> \endverbatim
56 *>
57 *> \param[in] QSIZ
58 *> \verbatim
59 *> QSIZ is INTEGER
60 *> The dimension of the orthogonal matrix used to reduce
61 *> the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
62 *> \endverbatim
63 *>
64 *> \param[in] N
65 *> \verbatim
66 *> N is INTEGER
67 *> The dimension of the symmetric tridiagonal matrix. N >= 0.
68 *> \endverbatim
69 *>
70 *> \param[in,out] D
71 *> \verbatim
72 *> D is DOUBLE PRECISION array, dimension (N)
73 *> On entry, the main diagonal of the tridiagonal matrix.
74 *> On exit, its eigenvalues.
75 *> \endverbatim
76 *>
77 *> \param[in] E
78 *> \verbatim
79 *> E is DOUBLE PRECISION array, dimension (N-1)
80 *> The off-diagonal elements of the tridiagonal matrix.
81 *> On exit, E has been destroyed.
82 *> \endverbatim
83 *>
84 *> \param[in,out] Q
85 *> \verbatim
86 *> Q is DOUBLE PRECISION array, dimension (LDQ, N)
87 *> On entry, Q must contain an N-by-N orthogonal matrix.
88 *> If ICOMPQ = 0 Q is not referenced.
89 *> If ICOMPQ = 1 On entry, Q is a subset of the columns of the
90 *> orthogonal matrix used to reduce the full
91 *> matrix to tridiagonal form corresponding to
92 *> the subset of the full matrix which is being
93 *> decomposed at this time.
94 *> If ICOMPQ = 2 On entry, Q will be the identity matrix.
95 *> On exit, Q contains the eigenvectors of the
96 *> tridiagonal matrix.
97 *> \endverbatim
98 *>
99 *> \param[in] LDQ
100 *> \verbatim
101 *> LDQ is INTEGER
102 *> The leading dimension of the array Q. If eigenvectors are
103 *> desired, then LDQ >= max(1,N). In any case, LDQ >= 1.
104 *> \endverbatim
105 *>
106 *> \param[out] QSTORE
107 *> \verbatim
108 *> QSTORE is DOUBLE PRECISION array, dimension (LDQS, N)
109 *> Referenced only when ICOMPQ = 1. Used to store parts of
110 *> the eigenvector matrix when the updating matrix multiplies
111 *> take place.
112 *> \endverbatim
113 *>
114 *> \param[in] LDQS
115 *> \verbatim
116 *> LDQS is INTEGER
117 *> The leading dimension of the array QSTORE. If ICOMPQ = 1,
118 *> then LDQS >= max(1,N). In any case, LDQS >= 1.
119 *> \endverbatim
120 *>
121 *> \param[out] WORK
122 *> \verbatim
123 *> WORK is DOUBLE PRECISION array,
124 *> If ICOMPQ = 0 or 1, the dimension of WORK must be at least
125 *> 1 + 3*N + 2*N*lg N + 3*N**2
126 *> ( lg( N ) = smallest integer k
127 *> such that 2^k >= N )
128 *> If ICOMPQ = 2, the dimension of WORK must be at least
129 *> 4*N + N**2.
130 *> \endverbatim
131 *>
132 *> \param[out] IWORK
133 *> \verbatim
134 *> IWORK is INTEGER array,
135 *> If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
136 *> 6 + 6*N + 5*N*lg N.
137 *> ( lg( N ) = smallest integer k
138 *> such that 2^k >= N )
139 *> If ICOMPQ = 2, the dimension of IWORK must be at least
140 *> 3 + 5*N.
141 *> \endverbatim
142 *>
143 *> \param[out] INFO
144 *> \verbatim
145 *> INFO is INTEGER
146 *> = 0: successful exit.
147 *> < 0: if INFO = -i, the i-th argument had an illegal value.
148 *> > 0: The algorithm failed to compute an eigenvalue while
149 *> working on the submatrix lying in rows and columns
150 *> INFO/(N+1) through mod(INFO,N+1).
151 *> \endverbatim
152 *
153 * Authors:
154 * ========
155 *
156 *> \author Univ. of Tennessee
157 *> \author Univ. of California Berkeley
158 *> \author Univ. of Colorado Denver
159 *> \author NAG Ltd.
160 *
161 *> \ingroup auxOTHERcomputational
162 *
163 *> \par Contributors:
164 * ==================
165 *>
166 *> Jeff Rutter, Computer Science Division, University of California
167 *> at Berkeley, USA
168 *
169 * =====================================================================
170  SUBROUTINE dlaed0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
171  $ WORK, IWORK, INFO )
172 *
173 * -- LAPACK computational routine --
174 * -- LAPACK is a software package provided by Univ. of Tennessee, --
175 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
176 *
177 * .. Scalar Arguments ..
178  INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
179 * ..
180 * .. Array Arguments ..
181  INTEGER IWORK( * )
182  DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
183  $ work( * )
184 * ..
185 *
186 * =====================================================================
187 *
188 * .. Parameters ..
189  DOUBLE PRECISION ZERO, ONE, TWO
190  parameter( zero = 0.d0, one = 1.d0, two = 2.d0 )
191 * ..
192 * .. Local Scalars ..
193  INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
194  $ igivpt, indxq, iperm, iprmpt, iq, iqptr, iwrem,
195  $ j, k, lgn, matsiz, msd2, smlsiz, smm1, spm1,
196  $ spm2, submat, subpbs, tlvls
197  DOUBLE PRECISION TEMP
198 * ..
199 * .. External Subroutines ..
200  EXTERNAL dcopy, dgemm, dlacpy, dlaed1, dlaed7, dsteqr,
201  $ xerbla
202 * ..
203 * .. External Functions ..
204  INTEGER ILAENV
205  EXTERNAL ilaenv
206 * ..
207 * .. Intrinsic Functions ..
208  INTRINSIC abs, dble, int, log, max
209 * ..
210 * .. Executable Statements ..
211 *
212 * Test the input parameters.
213 *
214  info = 0
215 *
216  IF( icompq.LT.0 .OR. icompq.GT.2 ) THEN
217  info = -1
218  ELSE IF( ( icompq.EQ.1 ) .AND. ( qsiz.LT.max( 0, n ) ) ) THEN
219  info = -2
220  ELSE IF( n.LT.0 ) THEN
221  info = -3
222  ELSE IF( ldq.LT.max( 1, n ) ) THEN
223  info = -7
224  ELSE IF( ldqs.LT.max( 1, n ) ) THEN
225  info = -9
226  END IF
227  IF( info.NE.0 ) THEN
228  CALL xerbla( 'DLAED0', -info )
229  RETURN
230  END IF
231 *
232 * Quick return if possible
233 *
234  IF( n.EQ.0 )
235  $ RETURN
236 *
237  smlsiz = ilaenv( 9, 'DLAED0', ' ', 0, 0, 0, 0 )
238 *
239 * Determine the size and placement of the submatrices, and save in
240 * the leading elements of IWORK.
241 *
242  iwork( 1 ) = n
243  subpbs = 1
244  tlvls = 0
245  10 CONTINUE
246  IF( iwork( subpbs ).GT.smlsiz ) THEN
247  DO 20 j = subpbs, 1, -1
248  iwork( 2*j ) = ( iwork( j )+1 ) / 2
249  iwork( 2*j-1 ) = iwork( j ) / 2
250  20 CONTINUE
251  tlvls = tlvls + 1
252  subpbs = 2*subpbs
253  GO TO 10
254  END IF
255  DO 30 j = 2, subpbs
256  iwork( j ) = iwork( j ) + iwork( j-1 )
257  30 CONTINUE
258 *
259 * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
260 * using rank-1 modifications (cuts).
261 *
262  spm1 = subpbs - 1
263  DO 40 i = 1, spm1
264  submat = iwork( i ) + 1
265  smm1 = submat - 1
266  d( smm1 ) = d( smm1 ) - abs( e( smm1 ) )
267  d( submat ) = d( submat ) - abs( e( smm1 ) )
268  40 CONTINUE
269 *
270  indxq = 4*n + 3
271  IF( icompq.NE.2 ) THEN
272 *
273 * Set up workspaces for eigenvalues only/accumulate new vectors
274 * routine
275 *
276  temp = log( dble( n ) ) / log( two )
277  lgn = int( temp )
278  IF( 2**lgn.LT.n )
279  $ lgn = lgn + 1
280  IF( 2**lgn.LT.n )
281  $ lgn = lgn + 1
282  iprmpt = indxq + n + 1
283  iperm = iprmpt + n*lgn
284  iqptr = iperm + n*lgn
285  igivpt = iqptr + n + 2
286  igivcl = igivpt + n*lgn
287 *
288  igivnm = 1
289  iq = igivnm + 2*n*lgn
290  iwrem = iq + n**2 + 1
291 *
292 * Initialize pointers
293 *
294  DO 50 i = 0, subpbs
295  iwork( iprmpt+i ) = 1
296  iwork( igivpt+i ) = 1
297  50 CONTINUE
298  iwork( iqptr ) = 1
299  END IF
300 *
301 * Solve each submatrix eigenproblem at the bottom of the divide and
302 * conquer tree.
303 *
304  curr = 0
305  DO 70 i = 0, spm1
306  IF( i.EQ.0 ) THEN
307  submat = 1
308  matsiz = iwork( 1 )
309  ELSE
310  submat = iwork( i ) + 1
311  matsiz = iwork( i+1 ) - iwork( i )
312  END IF
313  IF( icompq.EQ.2 ) THEN
314  CALL dsteqr( 'I', matsiz, d( submat ), e( submat ),
315  $ q( submat, submat ), ldq, work, info )
316  IF( info.NE.0 )
317  $ GO TO 130
318  ELSE
319  CALL dsteqr( 'I', matsiz, d( submat ), e( submat ),
320  $ work( iq-1+iwork( iqptr+curr ) ), matsiz, work,
321  $ info )
322  IF( info.NE.0 )
323  $ GO TO 130
324  IF( icompq.EQ.1 ) THEN
325  CALL dgemm( 'N', 'N', qsiz, matsiz, matsiz, one,
326  $ q( 1, submat ), ldq, work( iq-1+iwork( iqptr+
327  $ curr ) ), matsiz, zero, qstore( 1, submat ),
328  $ ldqs )
329  END IF
330  iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2
331  curr = curr + 1
332  END IF
333  k = 1
334  DO 60 j = submat, iwork( i+1 )
335  iwork( indxq+j ) = k
336  k = k + 1
337  60 CONTINUE
338  70 CONTINUE
339 *
340 * Successively merge eigensystems of adjacent submatrices
341 * into eigensystem for the corresponding larger matrix.
342 *
343 * while ( SUBPBS > 1 )
344 *
345  curlvl = 1
346  80 CONTINUE
347  IF( subpbs.GT.1 ) THEN
348  spm2 = subpbs - 2
349  DO 90 i = 0, spm2, 2
350  IF( i.EQ.0 ) THEN
351  submat = 1
352  matsiz = iwork( 2 )
353  msd2 = iwork( 1 )
354  curprb = 0
355  ELSE
356  submat = iwork( i ) + 1
357  matsiz = iwork( i+2 ) - iwork( i )
358  msd2 = matsiz / 2
359  curprb = curprb + 1
360  END IF
361 *
362 * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
363 * into an eigensystem of size MATSIZ.
364 * DLAED1 is used only for the full eigensystem of a tridiagonal
365 * matrix.
366 * DLAED7 handles the cases in which eigenvalues only or eigenvalues
367 * and eigenvectors of a full symmetric matrix (which was reduced to
368 * tridiagonal form) are desired.
369 *
370  IF( icompq.EQ.2 ) THEN
371  CALL dlaed1( matsiz, d( submat ), q( submat, submat ),
372  $ ldq, iwork( indxq+submat ),
373  $ e( submat+msd2-1 ), msd2, work,
374  $ iwork( subpbs+1 ), info )
375  ELSE
376  CALL dlaed7( icompq, matsiz, qsiz, tlvls, curlvl, curprb,
377  $ d( submat ), qstore( 1, submat ), ldqs,
378  $ iwork( indxq+submat ), e( submat+msd2-1 ),
379  $ msd2, work( iq ), iwork( iqptr ),
380  $ iwork( iprmpt ), iwork( iperm ),
381  $ iwork( igivpt ), iwork( igivcl ),
382  $ work( igivnm ), work( iwrem ),
383  $ iwork( subpbs+1 ), info )
384  END IF
385  IF( info.NE.0 )
386  $ GO TO 130
387  iwork( i / 2+1 ) = iwork( i+2 )
388  90 CONTINUE
389  subpbs = subpbs / 2
390  curlvl = curlvl + 1
391  GO TO 80
392  END IF
393 *
394 * end while
395 *
396 * Re-merge the eigenvalues/vectors which were deflated at the final
397 * merge step.
398 *
399  IF( icompq.EQ.1 ) THEN
400  DO 100 i = 1, n
401  j = iwork( indxq+i )
402  work( i ) = d( j )
403  CALL dcopy( qsiz, qstore( 1, j ), 1, q( 1, i ), 1 )
404  100 CONTINUE
405  CALL dcopy( n, work, 1, d, 1 )
406  ELSE IF( icompq.EQ.2 ) THEN
407  DO 110 i = 1, n
408  j = iwork( indxq+i )
409  work( i ) = d( j )
410  CALL dcopy( n, q( 1, j ), 1, work( n*i+1 ), 1 )
411  110 CONTINUE
412  CALL dcopy( n, work, 1, d, 1 )
413  CALL dlacpy( 'A', n, n, work( n+1 ), n, q, ldq )
414  ELSE
415  DO 120 i = 1, n
416  j = iwork( indxq+i )
417  work( i ) = d( j )
418  120 CONTINUE
419  CALL dcopy( n, work, 1, d, 1 )
420  END IF
421  GO TO 140
422 *
423  130 CONTINUE
424  info = submat*( n+1 ) + submat + matsiz - 1
425 *
426  140 CONTINUE
427  RETURN
428 *
429 * End of DLAED0
430 *
431  END
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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:131
subroutine dlaed1(N, D, Q, LDQ, INDXQ, RHO, CUTPNT, WORK, IWORK, INFO)
DLAED1 used by DSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a ...
Definition: dlaed1.f:163
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 dlaed7(ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q, LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK, INFO)
DLAED7 used by DSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a ...
Definition: dlaed7.f:260
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187