LAPACK: SRC/sstedc.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
00002      $                   LIWORK, INFO )
00003 *
00004 *  -- LAPACK driver routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          COMPZ
00011       INTEGER            INFO, LDZ, LIWORK, LWORK, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       INTEGER            IWORK( * )
00015       REAL               D( * ), E( * ), WORK( * ), Z( LDZ, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  SSTEDC computes all eigenvalues and, optionally, eigenvectors of a
00022 *  symmetric tridiagonal matrix using the divide and conquer method.
00023 *  The eigenvectors of a full or band real symmetric matrix can also be
00024 *  found if SSYTRD or SSPTRD or SSBTRD has been used to reduce this
00025 *  matrix to tridiagonal form.
00026 *
00027 *  This code makes very mild assumptions about floating point
00028 *  arithmetic. It will work on machines with a guard digit in
00029 *  add/subtract, or on those binary machines without guard digits
00030 *  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
00031 *  It could conceivably fail on hexadecimal or decimal machines
00032 *  without guard digits, but we know of none.  See SLAED3 for details.
00033 *
00034 *  Arguments
00035 *  =========
00036 *
00037 *  COMPZ   (input) CHARACTER*1
00038 *          = 'N':  Compute eigenvalues only.
00039 *          = 'I':  Compute eigenvectors of tridiagonal matrix also.
00040 *          = 'V':  Compute eigenvectors of original dense symmetric
00041 *                  matrix also.  On entry, Z contains the orthogonal
00042 *                  matrix used to reduce the original matrix to
00043 *                  tridiagonal form.
00044 *
00045 *  N       (input) INTEGER
00046 *          The dimension of the symmetric tridiagonal matrix.  N >= 0.
00047 *
00048 *  D       (input/output) REAL array, dimension (N)
00049 *          On entry, the diagonal elements of the tridiagonal matrix.
00050 *          On exit, if INFO = 0, the eigenvalues in ascending order.
00051 *
00052 *  E       (input/output) REAL array, dimension (N-1)
00053 *          On entry, the subdiagonal elements of the tridiagonal matrix.
00054 *          On exit, E has been destroyed.
00055 *
00056 *  Z       (input/output) REAL array, dimension (LDZ,N)
00057 *          On entry, if COMPZ = 'V', then Z contains the orthogonal
00058 *          matrix used in the reduction to tridiagonal form.
00059 *          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
00060 *          orthonormal eigenvectors of the original symmetric matrix,
00061 *          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
00062 *          of the symmetric tridiagonal matrix.
00063 *          If  COMPZ = 'N', then Z is not referenced.
00064 *
00065 *  LDZ     (input) INTEGER
00066 *          The leading dimension of the array Z.  LDZ >= 1.
00067 *          If eigenvectors are desired, then LDZ >= max(1,N).
00068 *
00069 *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
00070 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00071 *
00072 *  LWORK   (input) INTEGER
00073 *          The dimension of the array WORK.
00074 *          If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
00075 *          If COMPZ = 'V' and N > 1 then LWORK must be at least
00076 *                         ( 1 + 3*N + 2*N*lg N + 3*N**2 ),
00077 *                         where lg( N ) = smallest integer k such
00078 *                         that 2**k >= N.
00079 *          If COMPZ = 'I' and N > 1 then LWORK must be at least
00080 *                         ( 1 + 4*N + N**2 ).
00081 *          Note that for COMPZ = 'I' or 'V', then if N is less than or
00082 *          equal to the minimum divide size, usually 25, then LWORK need
00083 *          only be max(1,2*(N-1)).
00084 *
00085 *          If LWORK = -1, then a workspace query is assumed; the routine
00086 *          only calculates the optimal size of the WORK array, returns
00087 *          this value as the first entry of the WORK array, and no error
00088 *          message related to LWORK is issued by XERBLA.
00089 *
00090 *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
00091 *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
00092 *
00093 *  LIWORK  (input) INTEGER
00094 *          The dimension of the array IWORK.
00095 *          If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
00096 *          If COMPZ = 'V' and N > 1 then LIWORK must be at least
00097 *                         ( 6 + 6*N + 5*N*lg N ).
00098 *          If COMPZ = 'I' and N > 1 then LIWORK must be at least
00099 *                         ( 3 + 5*N ).
00100 *          Note that for COMPZ = 'I' or 'V', then if N is less than or
00101 *          equal to the minimum divide size, usually 25, then LIWORK
00102 *          need only be 1.
00103 *
00104 *          If LIWORK = -1, then a workspace query is assumed; the
00105 *          routine only calculates the optimal size of the IWORK array,
00106 *          returns this value as the first entry of the IWORK array, and
00107 *          no error message related to LIWORK is issued by XERBLA.
00108 *
00109 *  INFO    (output) INTEGER
00110 *          = 0:  successful exit.
00111 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00112 *          > 0:  The algorithm failed to compute an eigenvalue while
00113 *                working on the submatrix lying in rows and columns
00114 *                INFO/(N+1) through mod(INFO,N+1).
00115 *
00116 *  Further Details
00117 *  ===============
00118 *
00119 *  Based on contributions by
00120 *     Jeff Rutter, Computer Science Division, University of California
00121 *     at Berkeley, USA
00122 *  Modified by Francoise Tisseur, University of Tennessee.
00123 *
00124 *  =====================================================================
00125 *
00126 *     .. Parameters ..
00127       REAL               ZERO, ONE, TWO
00128       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
00129 *     ..
00130 *     .. Local Scalars ..
00131       LOGICAL            LQUERY
00132       INTEGER            FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
00133      $                   LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
00134       REAL               EPS, ORGNRM, P, TINY
00135 *     ..
00136 *     .. External Functions ..
00137       LOGICAL            LSAME
00138       INTEGER            ILAENV
00139       REAL               SLAMCH, SLANST
00140       EXTERNAL           ILAENV, LSAME, SLAMCH, SLANST
00141 *     ..
00142 *     .. External Subroutines ..
00143       EXTERNAL           SGEMM, SLACPY, SLAED0, SLASCL, SLASET, SLASRT,
00144      $                   SSTEQR, SSTERF, SSWAP, XERBLA
00145 *     ..
00146 *     .. Intrinsic Functions ..
00147       INTRINSIC          ABS, INT, LOG, MAX, MOD, REAL, SQRT
00148 *     ..
00149 *     .. Executable Statements ..
00150 *
00151 *     Test the input parameters.
00152 *
00153       INFO = 0
00154       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
00155 *
00156       IF( LSAME( COMPZ, 'N' ) ) THEN
00157          ICOMPZ = 0
00158       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
00159          ICOMPZ = 1
00160       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
00161          ICOMPZ = 2
00162       ELSE
00163          ICOMPZ = -1
00164       END IF
00165       IF( ICOMPZ.LT.0 ) THEN
00166          INFO = -1
00167       ELSE IF( N.LT.0 ) THEN
00168          INFO = -2
00169       ELSE IF( ( LDZ.LT.1 ) .OR.
00170      $         ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
00171          INFO = -6
00172       END IF
00173 *
00174       IF( INFO.EQ.0 ) THEN
00175 *
00176 *        Compute the workspace requirements
00177 *
00178          SMLSIZ = ILAENV( 9, 'SSTEDC', ' ', 0, 0, 0, 0 )
00179          IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN
00180             LIWMIN = 1
00181             LWMIN = 1
00182          ELSE IF( N.LE.SMLSIZ ) THEN
00183             LIWMIN = 1
00184             LWMIN = 2*( N - 1 )
00185          ELSE
00186             LGN = INT( LOG( REAL( N ) )/LOG( TWO ) )
00187             IF( 2**LGN.LT.N )
00188      $         LGN = LGN + 1
00189             IF( 2**LGN.LT.N )
00190      $         LGN = LGN + 1
00191             IF( ICOMPZ.EQ.1 ) THEN
00192                LWMIN = 1 + 3*N + 2*N*LGN + 3*N**2
00193                LIWMIN = 6 + 6*N + 5*N*LGN
00194             ELSE IF( ICOMPZ.EQ.2 ) THEN
00195                LWMIN = 1 + 4*N + N**2
00196                LIWMIN = 3 + 5*N
00197             END IF
00198          END IF
00199          WORK( 1 ) = LWMIN
00200          IWORK( 1 ) = LIWMIN
00201 *
00202          IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN
00203             INFO = -8
00204          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN
00205             INFO = -10
00206          END IF
00207       END IF
00208 *
00209       IF( INFO.NE.0 ) THEN
00210          CALL XERBLA( 'SSTEDC', -INFO )
00211          RETURN
00212       ELSE IF (LQUERY) THEN
00213          RETURN
00214       END IF
00215 *
00216 *     Quick return if possible
00217 *
00218       IF( N.EQ.0 )
00219      $   RETURN
00220       IF( N.EQ.1 ) THEN
00221          IF( ICOMPZ.NE.0 )
00222      $      Z( 1, 1 ) = ONE
00223          RETURN
00224       END IF
00225 *
00226 *     If the following conditional clause is removed, then the routine
00227 *     will use the Divide and Conquer routine to compute only the
00228 *     eigenvalues, which requires (3N + 3N**2) real workspace and
00229 *     (2 + 5N + 2N lg(N)) integer workspace.
00230 *     Since on many architectures SSTERF is much faster than any other
00231 *     algorithm for finding eigenvalues only, it is used here
00232 *     as the default. If the conditional clause is removed, then
00233 *     information on the size of workspace needs to be changed.
00234 *
00235 *     If COMPZ = 'N', use SSTERF to compute the eigenvalues.
00236 *
00237       IF( ICOMPZ.EQ.0 ) THEN
00238          CALL SSTERF( N, D, E, INFO )
00239          GO TO 50
00240       END IF
00241 *
00242 *     If N is smaller than the minimum divide size (SMLSIZ+1), then
00243 *     solve the problem with another solver.
00244 *
00245       IF( N.LE.SMLSIZ ) THEN
00246 *
00247          CALL SSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
00248 *
00249       ELSE
00250 *
00251 *        If COMPZ = 'V', the Z matrix must be stored elsewhere for later
00252 *        use.
00253 *
00254          IF( ICOMPZ.EQ.1 ) THEN
00255             STOREZ = 1 + N*N
00256          ELSE
00257             STOREZ = 1
00258          END IF
00259 *
00260          IF( ICOMPZ.EQ.2 ) THEN
00261             CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
00262          END IF
00263 *
00264 *        Scale.
00265 *
00266          ORGNRM = SLANST( 'M', N, D, E )
00267          IF( ORGNRM.EQ.ZERO )
00268      $      GO TO 50
00269 *
00270          EPS = SLAMCH( 'Epsilon' )
00271 *
00272          START = 1
00273 *
00274 *        while ( START <= N )
00275 *
00276    10    CONTINUE
00277          IF( START.LE.N ) THEN
00278 *
00279 *           Let FINISH be the position of the next subdiagonal entry
00280 *           such that E( FINISH ) <= TINY or FINISH = N if no such
00281 *           subdiagonal exists.  The matrix identified by the elements
00282 *           between START and FINISH constitutes an independent
00283 *           sub-problem.
00284 *
00285             FINISH = START
00286    20       CONTINUE
00287             IF( FINISH.LT.N ) THEN
00288                TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
00289      $                    SQRT( ABS( D( FINISH+1 ) ) )
00290                IF( ABS( E( FINISH ) ).GT.TINY ) THEN
00291                   FINISH = FINISH + 1
00292                   GO TO 20
00293                END IF
00294             END IF
00295 *
00296 *           (Sub) Problem determined.  Compute its size and solve it.
00297 *
00298             M = FINISH - START + 1
00299             IF( M.EQ.1 ) THEN
00300                START = FINISH + 1
00301                GO TO 10
00302             END IF
00303             IF( M.GT.SMLSIZ ) THEN
00304 *
00305 *              Scale.
00306 *
00307                ORGNRM = SLANST( 'M', M, D( START ), E( START ) )
00308                CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
00309      $                      INFO )
00310                CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
00311      $                      M-1, INFO )
00312 *
00313                IF( ICOMPZ.EQ.1 ) THEN
00314                   STRTRW = 1
00315                ELSE
00316                   STRTRW = START
00317                END IF
00318                CALL SLAED0( ICOMPZ, N, M, D( START ), E( START ),
00319      $                      Z( STRTRW, START ), LDZ, WORK( 1 ), N,
00320      $                      WORK( STOREZ ), IWORK, INFO )
00321                IF( INFO.NE.0 ) THEN
00322                   INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
00323      $                   MOD( INFO, ( M+1 ) ) + START - 1
00324                   GO TO 50
00325                END IF
00326 *
00327 *              Scale back.
00328 *
00329                CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
00330      $                      INFO )
00331 *
00332             ELSE
00333                IF( ICOMPZ.EQ.1 ) THEN
00334 *
00335 *                 Since QR won't update a Z matrix which is larger than
00336 *                 the length of D, we must solve the sub-problem in a
00337 *                 workspace and then multiply back into Z.
00338 *
00339                   CALL SSTEQR( 'I', M, D( START ), E( START ), WORK, M,
00340      $                         WORK( M*M+1 ), INFO )
00341                   CALL SLACPY( 'A', N, M, Z( 1, START ), LDZ,
00342      $                         WORK( STOREZ ), N )
00343                   CALL SGEMM( 'N', 'N', N, M, M, ONE,
00344      $                        WORK( STOREZ ), N, WORK, M, ZERO,
00345      $                        Z( 1, START ), LDZ )
00346                ELSE IF( ICOMPZ.EQ.2 ) THEN
00347                   CALL SSTEQR( 'I', M, D( START ), E( START ),
00348      $                         Z( START, START ), LDZ, WORK, INFO )
00349                ELSE
00350                   CALL SSTERF( M, D( START ), E( START ), INFO )
00351                END IF
00352                IF( INFO.NE.0 ) THEN
00353                   INFO = START*( N+1 ) + FINISH
00354                   GO TO 50
00355                END IF
00356             END IF
00357 *
00358             START = FINISH + 1
00359             GO TO 10
00360          END IF
00361 *
00362 *        endwhile
00363 *
00364 *        If the problem split any number of times, then the eigenvalues
00365 *        will not be properly ordered.  Here we permute the eigenvalues
00366 *        (and the associated eigenvectors) into ascending order.
00367 *
00368          IF( M.NE.N ) THEN
00369             IF( ICOMPZ.EQ.0 ) THEN
00370 *
00371 *              Use Quick Sort
00372 *
00373                CALL SLASRT( 'I', N, D, INFO )
00374 *
00375             ELSE
00376 *
00377 *              Use Selection Sort to minimize swaps of eigenvectors
00378 *
00379                DO 40 II = 2, N
00380                   I = II - 1
00381                   K = I
00382                   P = D( I )
00383                   DO 30 J = II, N
00384                      IF( D( J ).LT.P ) THEN
00385                         K = J
00386                         P = D( J )
00387                      END IF
00388    30             CONTINUE
00389                   IF( K.NE.I ) THEN
00390                      D( K ) = D( I )
00391                      D( I ) = P
00392                      CALL SSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
00393                   END IF
00394    40          CONTINUE
00395             END IF
00396          END IF
00397       END IF
00398 *
00399    50 CONTINUE
00400       WORK( 1 ) = LWMIN
00401       IWORK( 1 ) = LIWMIN
00402 *
00403       RETURN
00404 *
00405 *     End of SSTEDC
00406 *
00407       END