LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ sstedc()

subroutine sstedc ( character compz,
integer n,
real, dimension( * ) d,
real, dimension( * ) e,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

SSTEDC

Download SSTEDC + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SSTEDC computes all eigenvalues and, optionally, eigenvectors of a
!> symmetric tridiagonal matrix using the divide and conquer method.
!> The eigenvectors of a full or band real symmetric matrix can also be
!> found if SSYTRD or SSPTRD or SSBTRD has been used to reduce this
!> matrix to tridiagonal form.
!>
!>
Parameters
[in]COMPZ
!>          COMPZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only.
!>          = 'I':  Compute eigenvectors of tridiagonal matrix also.
!>          = 'V':  Compute eigenvectors of original dense symmetric
!>                  matrix also.  On entry, Z contains the orthogonal
!>                  matrix used to reduce the original matrix to
!>                  tridiagonal form.
!>
[in]N
!>          N is INTEGER
!>          The dimension of the symmetric tridiagonal matrix.  N >= 0.
!>
[in,out]D
!>          D is REAL array, dimension (N)
!>          On entry, the diagonal elements of the tridiagonal matrix.
!>          On exit, if INFO = 0, the eigenvalues in ascending order.
!>
[in,out]E
!>          E is REAL array, dimension (N-1)
!>          On entry, the subdiagonal elements of the tridiagonal matrix.
!>          On exit, E has been destroyed.
!>
[in,out]Z
!>          Z is REAL array, dimension (LDZ,N)
!>          On entry, if COMPZ = 'V', then Z contains the orthogonal
!>          matrix used in the reduction to tridiagonal form.
!>          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
!>          orthonormal eigenvectors of the original symmetric matrix,
!>          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
!>          of the symmetric tridiagonal matrix.
!>          If  COMPZ = 'N', then Z is not referenced.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1.
!>          If eigenvectors are desired, then LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is REAL array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
!>          If COMPZ = 'V' and N > 1 then LWORK must be at least
!>                         ( 1 + 3*N + 2*N*lg N + 4*N**2 ),
!>                         where lg( N ) = smallest integer k such
!>                         that 2**k >= N.
!>          If COMPZ = 'I' and N > 1 then LWORK must be at least
!>                         ( 1 + 4*N + N**2 ).
!>          Note that for COMPZ = 'I' or 'V', then if N is less than or
!>          equal to the minimum divide size, usually 25, then LWORK need
!>          only be max(1,2*(N-1)).
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.
!>          If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
!>          If COMPZ = 'V' and N > 1 then LIWORK must be at least
!>                         ( 6 + 6*N + 5*N*lg N ).
!>          If COMPZ = 'I' and N > 1 then LIWORK must be at least
!>                         ( 3 + 5*N ).
!>          Note that for COMPZ = 'I' or 'V', then if N is less than or
!>          equal to the minimum divide size, usually 25, then LIWORK
!>          need only be 1.
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal size of the IWORK array,
!>          returns this value as the first entry of the IWORK array, and
!>          no error message related to LIWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  The algorithm failed to compute an eigenvalue while
!>                working on the submatrix lying in rows and columns
!>                INFO/(N+1) through mod(INFO,N+1).
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

Definition at line 178 of file sstedc.f.

180*
181* -- LAPACK computational routine --
182* -- LAPACK is a software package provided by Univ. of Tennessee, --
183* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
184*
185* .. Scalar Arguments ..
186 CHARACTER COMPZ
187 INTEGER INFO, LDZ, LIWORK, LWORK, N
188* ..
189* .. Array Arguments ..
190 INTEGER IWORK( * )
191 REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
192* ..
193*
194* =====================================================================
195*
196* .. Parameters ..
197 REAL ZERO, ONE, TWO
198 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
199* ..
200* .. Local Scalars ..
201 LOGICAL LQUERY
202 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
203 $ LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
204 REAL EPS, ORGNRM, P, TINY
205* ..
206* .. External Functions ..
207 LOGICAL LSAME
208 INTEGER ILAENV
209 REAL SLAMCH, SLANST, SROUNDUP_LWORK
210 EXTERNAL ilaenv, lsame, slamch, slanst,
211 $ sroundup_lwork
212* ..
213* .. External Subroutines ..
214 EXTERNAL sgemm, slacpy, slaed0, slascl, slaset,
215 $ slasrt,
216 $ ssteqr, ssterf, sswap, xerbla
217* ..
218* .. Intrinsic Functions ..
219 INTRINSIC abs, int, log, max, mod, real, sqrt
220* ..
221* .. Executable Statements ..
222*
223* Test the input parameters.
224*
225 info = 0
226 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
227*
228 IF( lsame( compz, 'N' ) ) THEN
229 icompz = 0
230 ELSE IF( lsame( compz, 'V' ) ) THEN
231 icompz = 1
232 ELSE IF( lsame( compz, 'I' ) ) THEN
233 icompz = 2
234 ELSE
235 icompz = -1
236 END IF
237 IF( icompz.LT.0 ) THEN
238 info = -1
239 ELSE IF( n.LT.0 ) THEN
240 info = -2
241 ELSE IF( ( ldz.LT.1 ) .OR.
242 $ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) ) THEN
243 info = -6
244 END IF
245*
246 IF( info.EQ.0 ) THEN
247*
248* Compute the workspace requirements
249*
250 smlsiz = ilaenv( 9, 'SSTEDC', ' ', 0, 0, 0, 0 )
251 IF( n.LE.1 .OR. icompz.EQ.0 ) THEN
252 liwmin = 1
253 lwmin = 1
254 ELSE IF( n.LE.smlsiz ) THEN
255 liwmin = 1
256 lwmin = 2*( n - 1 )
257 ELSE
258 lgn = int( log( real( n ) )/log( two ) )
259 IF( 2**lgn.LT.n )
260 $ lgn = lgn + 1
261 IF( 2**lgn.LT.n )
262 $ lgn = lgn + 1
263 IF( icompz.EQ.1 ) THEN
264 lwmin = 1 + 3*n + 2*n*lgn + 4*n**2
265 liwmin = 6 + 6*n + 5*n*lgn
266 ELSE IF( icompz.EQ.2 ) THEN
267 lwmin = 1 + 4*n + n**2
268 liwmin = 3 + 5*n
269 END IF
270 END IF
271 work( 1 ) = sroundup_lwork(lwmin)
272 iwork( 1 ) = liwmin
273*
274 IF( lwork.LT.lwmin .AND. .NOT. lquery ) THEN
275 info = -8
276 ELSE IF( liwork.LT.liwmin .AND. .NOT. lquery ) THEN
277 info = -10
278 END IF
279 END IF
280*
281 IF( info.NE.0 ) THEN
282 CALL xerbla( 'SSTEDC', -info )
283 RETURN
284 ELSE IF (lquery) THEN
285 RETURN
286 END IF
287*
288* Quick return if possible
289*
290 IF( n.EQ.0 )
291 $ RETURN
292 IF( n.EQ.1 ) THEN
293 IF( icompz.NE.0 )
294 $ z( 1, 1 ) = one
295 RETURN
296 END IF
297*
298* If the following conditional clause is removed, then the routine
299* will use the Divide and Conquer routine to compute only the
300* eigenvalues, which requires (3N + 3N**2) real workspace and
301* (2 + 5N + 2N lg(N)) integer workspace.
302* Since on many architectures SSTERF is much faster than any other
303* algorithm for finding eigenvalues only, it is used here
304* as the default. If the conditional clause is removed, then
305* information on the size of workspace needs to be changed.
306*
307* If COMPZ = 'N', use SSTERF to compute the eigenvalues.
308*
309 IF( icompz.EQ.0 ) THEN
310 CALL ssterf( n, d, e, info )
311 GO TO 50
312 END IF
313*
314* If N is smaller than the minimum divide size (SMLSIZ+1), then
315* solve the problem with another solver.
316*
317 IF( n.LE.smlsiz ) THEN
318*
319 CALL ssteqr( compz, n, d, e, z, ldz, work, info )
320*
321 ELSE
322*
323* If COMPZ = 'V', the Z matrix must be stored elsewhere for later
324* use.
325*
326 IF( icompz.EQ.1 ) THEN
327 storez = 1 + n*n
328 ELSE
329 storez = 1
330 END IF
331*
332 IF( icompz.EQ.2 ) THEN
333 CALL slaset( 'Full', n, n, zero, one, z, ldz )
334 END IF
335*
336* Scale.
337*
338 orgnrm = slanst( 'M', n, d, e )
339 IF( orgnrm.EQ.zero )
340 $ GO TO 50
341*
342 eps = slamch( 'Epsilon' )
343*
344 start = 1
345*
346* while ( START <= N )
347*
348 10 CONTINUE
349 IF( start.LE.n ) THEN
350*
351* Let FINISH be the position of the next subdiagonal entry
352* such that E( FINISH ) <= TINY or FINISH = N if no such
353* subdiagonal exists. The matrix identified by the elements
354* between START and FINISH constitutes an independent
355* sub-problem.
356*
357 finish = start
358 20 CONTINUE
359 IF( finish.LT.n ) THEN
360 tiny = eps*sqrt( abs( d( finish ) ) )*
361 $ sqrt( abs( d( finish+1 ) ) )
362 IF( abs( e( finish ) ).GT.tiny ) THEN
363 finish = finish + 1
364 GO TO 20
365 END IF
366 END IF
367*
368* (Sub) Problem determined. Compute its size and solve it.
369*
370 m = finish - start + 1
371 IF( m.EQ.1 ) THEN
372 start = finish + 1
373 GO TO 10
374 END IF
375 IF( m.GT.smlsiz ) THEN
376*
377* Scale.
378*
379 orgnrm = slanst( 'M', m, d( start ), e( start ) )
380 CALL slascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ),
381 $ m,
382 $ info )
383 CALL slascl( 'G', 0, 0, orgnrm, one, m-1, 1,
384 $ e( start ),
385 $ m-1, info )
386*
387 IF( icompz.EQ.1 ) THEN
388 strtrw = 1
389 ELSE
390 strtrw = start
391 END IF
392 CALL slaed0( icompz, n, m, d( start ), e( start ),
393 $ z( strtrw, start ), ldz, work( 1 ), n,
394 $ work( storez ), iwork, info )
395 IF( info.NE.0 ) THEN
396 info = ( info / ( m+1 )+start-1 )*( n+1 ) +
397 $ mod( info, ( m+1 ) ) + start - 1
398 GO TO 50
399 END IF
400*
401* Scale back.
402*
403 CALL slascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ),
404 $ m,
405 $ info )
406*
407 ELSE
408 IF( icompz.EQ.1 ) THEN
409*
410* Since QR won't update a Z matrix which is larger than
411* the length of D, we must solve the sub-problem in a
412* workspace and then multiply back into Z.
413*
414 CALL ssteqr( 'I', m, d( start ), e( start ), work,
415 $ m,
416 $ work( m*m+1 ), info )
417 CALL slacpy( 'A', n, m, z( 1, start ), ldz,
418 $ work( storez ), n )
419 CALL sgemm( 'N', 'N', n, m, m, one,
420 $ work( storez ), n, work, m, zero,
421 $ z( 1, start ), ldz )
422 ELSE IF( icompz.EQ.2 ) THEN
423 CALL ssteqr( 'I', m, d( start ), e( start ),
424 $ z( start, start ), ldz, work, info )
425 ELSE
426 CALL ssterf( m, d( start ), e( start ), info )
427 END IF
428 IF( info.NE.0 ) THEN
429 info = start*( n+1 ) + finish
430 GO TO 50
431 END IF
432 END IF
433*
434 start = finish + 1
435 GO TO 10
436 END IF
437*
438* endwhile
439*
440 IF( icompz.EQ.0 ) THEN
441*
442* Use Quick Sort
443*
444 CALL slasrt( 'I', n, d, info )
445*
446 ELSE
447*
448* Use Selection Sort to minimize swaps of eigenvectors
449*
450 DO 40 ii = 2, n
451 i = ii - 1
452 k = i
453 p = d( i )
454 DO 30 j = ii, n
455 IF( d( j ).LT.p ) THEN
456 k = j
457 p = d( j )
458 END IF
459 30 CONTINUE
460 IF( k.NE.i ) THEN
461 d( k ) = d( i )
462 d( i ) = p
463 CALL sswap( n, z( 1, i ), 1, z( 1, k ), 1 )
464 END IF
465 40 CONTINUE
466 END IF
467 END IF
468*
469 50 CONTINUE
470 work( 1 ) = sroundup_lwork(lwmin)
471 iwork( 1 ) = liwmin
472*
473 RETURN
474*
475* End of SSTEDC
476*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:101
slaed0
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:170
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slanst
real function slanst(norm, n, d, e)
SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slanst.f:98
slascl
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:142
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:108
slasrt
subroutine slasrt(id, n, d, info)
SLASRT sorts numbers in increasing or decreasing order.
Definition slasrt.f:86
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sroundup_lwork
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
Definition sroundup_lwork.f:59
ssteqr
subroutine ssteqr(compz, n, d, e, z, ldz, work, info)
SSTEQR
Definition ssteqr.f:129
ssterf
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:84
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: