LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine slaed1(N, D, Q, LDQ, INDXQ, RHO, CUTPNT, WORK, IWORK, INFO)
SLAED1 used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a ...
Definition: slaed1.f:165
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:133
subroutine slaed0(ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, WORK, IWORK, INFO)
SLAED0 used by sstedc. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmet...
Definition: slaed0.f:174
subroutine slaed7(ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q, LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK, INFO)
SLAED7 used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a ...
Definition: slaed7.f:262
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53