LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
slaed0.f
Go to the documentation of this file.
1 *> \brief \b SLAED0 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 SLAED0 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaed0.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaed0.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaed0.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLAED0( 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 * REAL D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
30 * $ WORK( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SLAED0 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 *> \date September 2012
162 *
163 *> \ingroup auxOTHERcomputational
164 *
165 *> \par Contributors:
166 * ==================
167 *>
168 *> Jeff Rutter, Computer Science Division, University of California
169 *> at Berkeley, USA
170 *
171 * =====================================================================
172  SUBROUTINE slaed0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
173  $ work, iwork, info )
174 *
175 * -- LAPACK computational routine (version 3.4.2) --
176 * -- LAPACK is a software package provided by Univ. of Tennessee, --
177 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
178 * September 2012
179 *
180 * .. Scalar Arguments ..
181  INTEGER icompq, info, ldq, ldqs, n, qsiz
182 * ..
183 * .. Array Arguments ..
184  INTEGER iwork( * )
185  REAL d( * ), e( * ), q( ldq, * ), qstore( ldqs, * ),
186  $ work( * )
187 * ..
188 *
189 * =====================================================================
190 *
191 * .. Parameters ..
192  REAL zero, one, two
193  parameter( zero = 0.e0, one = 1.e0, two = 2.e0 )
194 * ..
195 * .. Local Scalars ..
196  INTEGER curlvl, curprb, curr, i, igivcl, igivnm,
197  $ igivpt, indxq, iperm, iprmpt, iq, iqptr, iwrem,
198  $ j, k, lgn, matsiz, msd2, smlsiz, smm1, spm1,
199  $ spm2, submat, subpbs, tlvls
200  REAL temp
201 * ..
202 * .. External Subroutines ..
203  EXTERNAL scopy, sgemm, slacpy, slaed1, slaed7, ssteqr,
204  $ xerbla
205 * ..
206 * .. External Functions ..
207  INTEGER ilaenv
208  EXTERNAL ilaenv
209 * ..
210 * .. Intrinsic Functions ..
211  INTRINSIC abs, int, log, max, real
212 * ..
213 * .. Executable Statements ..
214 *
215 * Test the input parameters.
216 *
217  info = 0
218 *
219  IF( icompq.LT.0 .OR. icompq.GT.2 ) THEN
220  info = -1
221  ELSE IF( ( icompq.EQ.1 ) .AND. ( qsiz.LT.max( 0, n ) ) ) THEN
222  info = -2
223  ELSE IF( n.LT.0 ) THEN
224  info = -3
225  ELSE IF( ldq.LT.max( 1, n ) ) THEN
226  info = -7
227  ELSE IF( ldqs.LT.max( 1, n ) ) THEN
228  info = -9
229  END IF
230  IF( info.NE.0 ) THEN
231  CALL xerbla( 'SLAED0', -info )
232  return
233  END IF
234 *
235 * Quick return if possible
236 *
237  IF( n.EQ.0 )
238  $ return
239 *
240  smlsiz = ilaenv( 9, 'SLAED0', ' ', 0, 0, 0, 0 )
241 *
242 * Determine the size and placement of the submatrices, and save in
243 * the leading elements of IWORK.
244 *
245  iwork( 1 ) = n
246  subpbs = 1
247  tlvls = 0
248  10 continue
249  IF( iwork( subpbs ).GT.smlsiz ) THEN
250  DO 20 j = subpbs, 1, -1
251  iwork( 2*j ) = ( iwork( j )+1 ) / 2
252  iwork( 2*j-1 ) = iwork( j ) / 2
253  20 continue
254  tlvls = tlvls + 1
255  subpbs = 2*subpbs
256  go to 10
257  END IF
258  DO 30 j = 2, subpbs
259  iwork( j ) = iwork( j ) + iwork( j-1 )
260  30 continue
261 *
262 * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
263 * using rank-1 modifications (cuts).
264 *
265  spm1 = subpbs - 1
266  DO 40 i = 1, spm1
267  submat = iwork( i ) + 1
268  smm1 = submat - 1
269  d( smm1 ) = d( smm1 ) - abs( e( smm1 ) )
270  d( submat ) = d( submat ) - abs( e( smm1 ) )
271  40 continue
272 *
273  indxq = 4*n + 3
274  IF( icompq.NE.2 ) THEN
275 *
276 * Set up workspaces for eigenvalues only/accumulate new vectors
277 * routine
278 *
279  temp = log( REAL( N ) ) / log( two )
280  lgn = int( temp )
281  IF( 2**lgn.LT.n )
282  $ lgn = lgn + 1
283  IF( 2**lgn.LT.n )
284  $ lgn = lgn + 1
285  iprmpt = indxq + n + 1
286  iperm = iprmpt + n*lgn
287  iqptr = iperm + n*lgn
288  igivpt = iqptr + n + 2
289  igivcl = igivpt + n*lgn
290 *
291  igivnm = 1
292  iq = igivnm + 2*n*lgn
293  iwrem = iq + n**2 + 1
294 *
295 * Initialize pointers
296 *
297  DO 50 i = 0, subpbs
298  iwork( iprmpt+i ) = 1
299  iwork( igivpt+i ) = 1
300  50 continue
301  iwork( iqptr ) = 1
302  END IF
303 *
304 * Solve each submatrix eigenproblem at the bottom of the divide and
305 * conquer tree.
306 *
307  curr = 0
308  DO 70 i = 0, spm1
309  IF( i.EQ.0 ) THEN
310  submat = 1
311  matsiz = iwork( 1 )
312  ELSE
313  submat = iwork( i ) + 1
314  matsiz = iwork( i+1 ) - iwork( i )
315  END IF
316  IF( icompq.EQ.2 ) THEN
317  CALL ssteqr( 'I', matsiz, d( submat ), e( submat ),
318  $ q( submat, submat ), ldq, work, info )
319  IF( info.NE.0 )
320  $ go to 130
321  ELSE
322  CALL ssteqr( 'I', matsiz, d( submat ), e( submat ),
323  $ work( iq-1+iwork( iqptr+curr ) ), matsiz, work,
324  $ info )
325  IF( info.NE.0 )
326  $ go to 130
327  IF( icompq.EQ.1 ) THEN
328  CALL sgemm( 'N', 'N', qsiz, matsiz, matsiz, one,
329  $ q( 1, submat ), ldq, work( iq-1+iwork( iqptr+
330  $ curr ) ), matsiz, zero, qstore( 1, submat ),
331  $ ldqs )
332  END IF
333  iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2
334  curr = curr + 1
335  END IF
336  k = 1
337  DO 60 j = submat, iwork( i+1 )
338  iwork( indxq+j ) = k
339  k = k + 1
340  60 continue
341  70 continue
342 *
343 * Successively merge eigensystems of adjacent submatrices
344 * into eigensystem for the corresponding larger matrix.
345 *
346 * while ( SUBPBS > 1 )
347 *
348  curlvl = 1
349  80 continue
350  IF( subpbs.GT.1 ) THEN
351  spm2 = subpbs - 2
352  DO 90 i = 0, spm2, 2
353  IF( i.EQ.0 ) THEN
354  submat = 1
355  matsiz = iwork( 2 )
356  msd2 = iwork( 1 )
357  curprb = 0
358  ELSE
359  submat = iwork( i ) + 1
360  matsiz = iwork( i+2 ) - iwork( i )
361  msd2 = matsiz / 2
362  curprb = curprb + 1
363  END IF
364 *
365 * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
366 * into an eigensystem of size MATSIZ.
367 * SLAED1 is used only for the full eigensystem of a tridiagonal
368 * matrix.
369 * SLAED7 handles the cases in which eigenvalues only or eigenvalues
370 * and eigenvectors of a full symmetric matrix (which was reduced to
371 * tridiagonal form) are desired.
372 *
373  IF( icompq.EQ.2 ) THEN
374  CALL slaed1( matsiz, d( submat ), q( submat, submat ),
375  $ ldq, iwork( indxq+submat ),
376  $ e( submat+msd2-1 ), msd2, work,
377  $ iwork( subpbs+1 ), info )
378  ELSE
379  CALL slaed7( icompq, matsiz, qsiz, tlvls, curlvl, curprb,
380  $ d( submat ), qstore( 1, submat ), ldqs,
381  $ iwork( indxq+submat ), e( submat+msd2-1 ),
382  $ msd2, work( iq ), iwork( iqptr ),
383  $ iwork( iprmpt ), iwork( iperm ),
384  $ iwork( igivpt ), iwork( igivcl ),
385  $ work( igivnm ), work( iwrem ),
386  $ iwork( subpbs+1 ), info )
387  END IF
388  IF( info.NE.0 )
389  $ go to 130
390  iwork( i / 2+1 ) = iwork( i+2 )
391  90 continue
392  subpbs = subpbs / 2
393  curlvl = curlvl + 1
394  go to 80
395  END IF
396 *
397 * end while
398 *
399 * Re-merge the eigenvalues/vectors which were deflated at the final
400 * merge step.
401 *
402  IF( icompq.EQ.1 ) THEN
403  DO 100 i = 1, n
404  j = iwork( indxq+i )
405  work( i ) = d( j )
406  CALL scopy( qsiz, qstore( 1, j ), 1, q( 1, i ), 1 )
407  100 continue
408  CALL scopy( n, work, 1, d, 1 )
409  ELSE IF( icompq.EQ.2 ) THEN
410  DO 110 i = 1, n
411  j = iwork( indxq+i )
412  work( i ) = d( j )
413  CALL scopy( n, q( 1, j ), 1, work( n*i+1 ), 1 )
414  110 continue
415  CALL scopy( n, work, 1, d, 1 )
416  CALL slacpy( 'A', n, n, work( n+1 ), n, q, ldq )
417  ELSE
418  DO 120 i = 1, n
419  j = iwork( indxq+i )
420  work( i ) = d( j )
421  120 continue
422  CALL scopy( n, work, 1, d, 1 )
423  END IF
424  go to 140
425 *
426  130 continue
427  info = submat*( n+1 ) + submat + matsiz - 1
428 *
429  140 continue
430  return
431 *
432 * End of SLAED0
433 *
434  END