LAPACK 3.3.1
Linear Algebra PACKage

ssyevd.f

Go to the documentation of this file.
00001       SUBROUTINE SSYEVD( JOBZ, UPLO, N, A, LDA, W, 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          JOBZ, UPLO
00011       INTEGER            INFO, LDA, LIWORK, LWORK, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       INTEGER            IWORK( * )
00015       REAL               A( LDA, * ), W( * ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  SSYEVD computes all eigenvalues and, optionally, eigenvectors of a
00022 *  real symmetric matrix A. If eigenvectors are desired, it uses a
00023 *  divide and conquer algorithm.
00024 *
00025 *  The divide and conquer algorithm makes very mild assumptions about
00026 *  floating point arithmetic. It will work on machines with a guard
00027 *  digit in add/subtract, or on those binary machines without guard
00028 *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
00029 *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
00030 *  without guard digits, but we know of none.
00031 *
00032 *  Because of large use of BLAS of level 3, SSYEVD needs N**2 more
00033 *  workspace than SSYEVX.
00034 *
00035 *  Arguments
00036 *  =========
00037 *
00038 *  JOBZ    (input) CHARACTER*1
00039 *          = 'N':  Compute eigenvalues only;
00040 *          = 'V':  Compute eigenvalues and eigenvectors.
00041 *
00042 *  UPLO    (input) CHARACTER*1
00043 *          = 'U':  Upper triangle of A is stored;
00044 *          = 'L':  Lower triangle of A is stored.
00045 *
00046 *  N       (input) INTEGER
00047 *          The order of the matrix A.  N >= 0.
00048 *
00049 *  A       (input/output) REAL array, dimension (LDA, N)
00050 *          On entry, the symmetric matrix A.  If UPLO = 'U', the
00051 *          leading N-by-N upper triangular part of A contains the
00052 *          upper triangular part of the matrix A.  If UPLO = 'L',
00053 *          the leading N-by-N lower triangular part of A contains
00054 *          the lower triangular part of the matrix A.
00055 *          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
00056 *          orthonormal eigenvectors of the matrix A.
00057 *          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
00058 *          or the upper triangle (if UPLO='U') of A, including the
00059 *          diagonal, is destroyed.
00060 *
00061 *  LDA     (input) INTEGER
00062 *          The leading dimension of the array A.  LDA >= max(1,N).
00063 *
00064 *  W       (output) REAL array, dimension (N)
00065 *          If INFO = 0, the eigenvalues in ascending order.
00066 *
00067 *  WORK    (workspace/output) REAL array,
00068 *                                         dimension (LWORK)
00069 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00070 *
00071 *  LWORK   (input) INTEGER
00072 *          The dimension of the array WORK.
00073 *          If N <= 1,               LWORK must be at least 1.
00074 *          If JOBZ = 'N' and N > 1, LWORK must be at least 2*N+1.
00075 *          If JOBZ = 'V' and N > 1, LWORK must be at least 
00076 *                                                1 + 6*N + 2*N**2.
00077 *
00078 *          If LWORK = -1, then a workspace query is assumed; the routine
00079 *          only calculates the optimal sizes of the WORK and IWORK
00080 *          arrays, returns these values as the first entries of the WORK
00081 *          and IWORK arrays, and no error message related to LWORK or
00082 *          LIWORK is issued by XERBLA.
00083 *
00084 *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
00085 *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
00086 *
00087 *  LIWORK  (input) INTEGER
00088 *          The dimension of the array IWORK.
00089 *          If N <= 1,                LIWORK must be at least 1.
00090 *          If JOBZ  = 'N' and N > 1, LIWORK must be at least 1.
00091 *          If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
00092 *
00093 *          If LIWORK = -1, then a workspace query is assumed; the
00094 *          routine only calculates the optimal sizes of the WORK and
00095 *          IWORK arrays, returns these values as the first entries of
00096 *          the WORK and IWORK arrays, and no error message related to
00097 *          LWORK or LIWORK is issued by XERBLA.
00098 *
00099 *  INFO    (output) INTEGER
00100 *          = 0:  successful exit
00101 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00102 *          > 0:  if INFO = i and JOBZ = 'N', then the algorithm failed
00103 *                to converge; i off-diagonal elements of an intermediate
00104 *                tridiagonal form did not converge to zero;
00105 *                if INFO = i and JOBZ = 'V', then the algorithm failed
00106 *                to compute an eigenvalue while working on the submatrix
00107 *                lying in rows and columns INFO/(N+1) through
00108 *                mod(INFO,N+1).
00109 *
00110 *  Further Details
00111 *  ===============
00112 *
00113 *  Based on contributions by
00114 *     Jeff Rutter, Computer Science Division, University of California
00115 *     at Berkeley, USA
00116 *  Modified by Francoise Tisseur, University of Tennessee.
00117 *
00118 *  Modified description of INFO. Sven, 16 Feb 05.
00119 *  =====================================================================
00120 *
00121 *     .. Parameters ..
00122       REAL               ZERO, ONE
00123       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00124 *     ..
00125 *     .. Local Scalars ..
00126 *
00127       LOGICAL            LOWER, LQUERY, WANTZ
00128       INTEGER            IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE,
00129      $                   LIOPT, LIWMIN, LLWORK, LLWRK2, LOPT, LWMIN
00130       REAL               ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
00131      $                   SMLNUM
00132 *     ..
00133 *     .. External Functions ..
00134       LOGICAL            LSAME
00135       INTEGER            ILAENV
00136       REAL               SLAMCH, SLANSY
00137       EXTERNAL           ILAENV, LSAME, SLAMCH, SLANSY
00138 *     ..
00139 *     .. External Subroutines ..
00140       EXTERNAL           SLACPY, SLASCL, SORMTR, SSCAL, SSTEDC, SSTERF,
00141      $                   SSYTRD, XERBLA
00142 *     ..
00143 *     .. Intrinsic Functions ..
00144       INTRINSIC          MAX, SQRT
00145 *     ..
00146 *     .. Executable Statements ..
00147 *
00148 *     Test the input parameters.
00149 *
00150       WANTZ = LSAME( JOBZ, 'V' )
00151       LOWER = LSAME( UPLO, 'L' )
00152       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
00153 *
00154       INFO = 0
00155       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
00156          INFO = -1
00157       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
00158          INFO = -2
00159       ELSE IF( N.LT.0 ) THEN
00160          INFO = -3
00161       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00162          INFO = -5
00163       END IF
00164 *
00165       IF( INFO.EQ.0 ) THEN
00166          IF( N.LE.1 ) THEN
00167             LIWMIN = 1
00168             LWMIN = 1
00169             LOPT = LWMIN
00170             LIOPT = LIWMIN
00171          ELSE
00172             IF( WANTZ ) THEN
00173                LIWMIN = 3 + 5*N
00174                LWMIN = 1 + 6*N + 2*N**2
00175             ELSE
00176                LIWMIN = 1
00177                LWMIN = 2*N + 1
00178             END IF
00179             LOPT = MAX( LWMIN, 2*N +
00180      $                  ILAENV( 1, 'SSYTRD', UPLO, N, -1, -1, -1 ) )
00181             LIOPT = LIWMIN
00182          END IF
00183          WORK( 1 ) = LOPT
00184          IWORK( 1 ) = LIOPT
00185 *
00186          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00187             INFO = -8
00188          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
00189             INFO = -10
00190          END IF
00191       END IF
00192 *
00193       IF( INFO.NE.0 ) THEN
00194          CALL XERBLA( 'SSYEVD', -INFO )
00195          RETURN
00196       ELSE IF( LQUERY ) THEN
00197          RETURN 
00198       END IF
00199 *
00200 *     Quick return if possible
00201 *
00202       IF( N.EQ.0 )
00203      $   RETURN
00204 *
00205       IF( N.EQ.1 ) THEN
00206          W( 1 ) = A( 1, 1 )
00207          IF( WANTZ )
00208      $      A( 1, 1 ) = ONE
00209          RETURN 
00210       END IF
00211 *
00212 *     Get machine constants.
00213 *
00214       SAFMIN = SLAMCH( 'Safe minimum' )
00215       EPS = SLAMCH( 'Precision' )
00216       SMLNUM = SAFMIN / EPS
00217       BIGNUM = ONE / SMLNUM
00218       RMIN = SQRT( SMLNUM )
00219       RMAX = SQRT( BIGNUM )
00220 *
00221 *     Scale matrix to allowable range, if necessary.
00222 *
00223       ANRM = SLANSY( 'M', UPLO, N, A, LDA, WORK )
00224       ISCALE = 0
00225       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
00226          ISCALE = 1
00227          SIGMA = RMIN / ANRM
00228       ELSE IF( ANRM.GT.RMAX ) THEN
00229          ISCALE = 1
00230          SIGMA = RMAX / ANRM
00231       END IF
00232       IF( ISCALE.EQ.1 )
00233      $   CALL SLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
00234 *
00235 *     Call SSYTRD to reduce symmetric matrix to tridiagonal form.
00236 *
00237       INDE = 1
00238       INDTAU = INDE + N
00239       INDWRK = INDTAU + N
00240       LLWORK = LWORK - INDWRK + 1
00241       INDWK2 = INDWRK + N*N
00242       LLWRK2 = LWORK - INDWK2 + 1
00243 *
00244       CALL SSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),
00245      $             WORK( INDWRK ), LLWORK, IINFO )
00246 *
00247 *     For eigenvalues only, call SSTERF.  For eigenvectors, first call
00248 *     SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
00249 *     tridiagonal matrix, then call SORMTR to multiply it by the
00250 *     Householder transformations stored in A.
00251 *
00252       IF( .NOT.WANTZ ) THEN
00253          CALL SSTERF( N, W, WORK( INDE ), INFO )
00254       ELSE
00255          CALL SSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
00256      $                WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
00257          CALL SORMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ),
00258      $                WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO )
00259          CALL SLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA )
00260       END IF
00261 *
00262 *     If matrix was scaled, then rescale eigenvalues appropriately.
00263 *
00264       IF( ISCALE.EQ.1 )
00265      $   CALL SSCAL( N, ONE / SIGMA, W, 1 )
00266 *
00267       WORK( 1 ) = LOPT
00268       IWORK( 1 ) = LIOPT
00269 *
00270       RETURN
00271 *
00272 *     End of SSYEVD
00273 *
00274       END
 All Files Functions