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