LAPACK 3.12.0
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 180 of file sstedc.f.

182*
183* -- LAPACK computational routine --
184* -- LAPACK is a software package provided by Univ. of Tennessee, --
185* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
186*
187* .. Scalar Arguments ..
188 CHARACTER COMPZ
189 INTEGER INFO, LDZ, LIWORK, LWORK, N
190* ..
191* .. Array Arguments ..
192 INTEGER IWORK( * )
193 REAL D( * ), E( * ), WORK( * ), Z( LDZ, * )
194* ..
195*
196* =====================================================================
197*
198* .. Parameters ..
199 REAL ZERO, ONE, TWO
200 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
201* ..
202* .. Local Scalars ..
203 LOGICAL LQUERY
204 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
205 $ LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
206 REAL EPS, ORGNRM, P, TINY
207* ..
208* .. External Functions ..
209 LOGICAL LSAME
210 INTEGER ILAENV
211 REAL SLAMCH, SLANST, SROUNDUP_LWORK
213* ..
214* .. External Subroutines ..
215 EXTERNAL sgemm, slacpy, slaed0, slascl, slaset, slasrt,
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 ), m,
381 $ info )
382 CALL slascl( 'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
383 $ m-1, info )
384*
385 IF( icompz.EQ.1 ) THEN
386 strtrw = 1
387 ELSE
388 strtrw = start
389 END IF
390 CALL slaed0( icompz, n, m, d( start ), e( start ),
391 $ z( strtrw, start ), ldz, work( 1 ), n,
392 $ work( storez ), iwork, info )
393 IF( info.NE.0 ) THEN
394 info = ( info / ( m+1 )+start-1 )*( n+1 ) +
395 $ mod( info, ( m+1 ) ) + start - 1
396 GO TO 50
397 END IF
398*
399* Scale back.
400*
401 CALL slascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
402 $ info )
403*
404 ELSE
405 IF( icompz.EQ.1 ) THEN
406*
407* Since QR won't update a Z matrix which is larger than
408* the length of D, we must solve the sub-problem in a
409* workspace and then multiply back into Z.
410*
411 CALL ssteqr( 'I', m, d( start ), e( start ), work, m,
412 $ work( m*m+1 ), info )
413 CALL slacpy( 'A', n, m, z( 1, start ), ldz,
414 $ work( storez ), n )
415 CALL sgemm( 'N', 'N', n, m, m, one,
416 $ work( storez ), n, work, m, zero,
417 $ z( 1, start ), ldz )
418 ELSE IF( icompz.EQ.2 ) THEN
419 CALL ssteqr( 'I', m, d( start ), e( start ),
420 $ z( start, start ), ldz, work, info )
421 ELSE
422 CALL ssterf( m, d( start ), e( start ), info )
423 END IF
424 IF( info.NE.0 ) THEN
425 info = start*( n+1 ) + finish
426 GO TO 50
427 END IF
428 END IF
429*
430 start = finish + 1
431 GO TO 10
432 END IF
433*
434* endwhile
435*
436 IF( icompz.EQ.0 ) THEN
437*
438* Use Quick Sort
439*
440 CALL slasrt( 'I', n, d, info )
441*
442 ELSE
443*
444* Use Selection Sort to minimize swaps of eigenvectors
445*
446 DO 40 ii = 2, n
447 i = ii - 1
448 k = i
449 p = d( i )
450 DO 30 j = ii, n
451 IF( d( j ).LT.p ) THEN
452 k = j
453 p = d( j )
454 END IF
455 30 CONTINUE
456 IF( k.NE.i ) THEN
457 d( k ) = d( i )
458 d( i ) = p
459 CALL sswap( n, z( 1, i ), 1, z( 1, k ), 1 )
460 END IF
461 40 CONTINUE
462 END IF
463 END IF
464*
465 50 CONTINUE
466 work( 1 ) = sroundup_lwork(lwmin)
467 iwork( 1 ) = liwmin
468*
469 RETURN
470*
471* End of SSTEDC
472*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
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:172
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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:100
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:143
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:110
subroutine slasrt(id, n, d, info)
SLASRT sorts numbers in increasing or decreasing order.
Definition slasrt.f:88
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine ssteqr(compz, n, d, e, z, ldz, work, info)
SSTEQR
Definition ssteqr.f:131
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:86
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: