LAPACK: SRC/dspev.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSPEV( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, INFO )
00002 *
00003 *  -- LAPACK driver routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          JOBZ, UPLO
00010       INTEGER            INFO, LDZ, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   AP( * ), W( * ), WORK( * ), Z( LDZ, * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DSPEV computes all the eigenvalues and, optionally, eigenvectors of a
00020 *  real symmetric matrix A in packed storage.
00021 *
00022 *  Arguments
00023 *  =========
00024 *
00025 *  JOBZ    (input) CHARACTER*1
00026 *          = 'N':  Compute eigenvalues only;
00027 *          = 'V':  Compute eigenvalues and eigenvectors.
00028 *
00029 *  UPLO    (input) CHARACTER*1
00030 *          = 'U':  Upper triangle of A is stored;
00031 *          = 'L':  Lower triangle of A is stored.
00032 *
00033 *  N       (input) INTEGER
00034 *          The order of the matrix A.  N >= 0.
00035 *
00036 *  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
00037 *          On entry, the upper or lower triangle of the symmetric matrix
00038 *          A, packed columnwise in a linear array.  The j-th column of A
00039 *          is stored in the array AP as follows:
00040 *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
00041 *          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
00042 *
00043 *          On exit, AP is overwritten by values generated during the
00044 *          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
00045 *          and first superdiagonal of the tridiagonal matrix T overwrite
00046 *          the corresponding elements of A, and if UPLO = 'L', the
00047 *          diagonal and first subdiagonal of T overwrite the
00048 *          corresponding elements of A.
00049 *
00050 *  W       (output) DOUBLE PRECISION array, dimension (N)
00051 *          If INFO = 0, the eigenvalues in ascending order.
00052 *
00053 *  Z       (output) DOUBLE PRECISION array, dimension (LDZ, N)
00054 *          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
00055 *          eigenvectors of the matrix A, with the i-th column of Z
00056 *          holding the eigenvector associated with W(i).
00057 *          If JOBZ = 'N', then Z is not referenced.
00058 *
00059 *  LDZ     (input) INTEGER
00060 *          The leading dimension of the array Z.  LDZ >= 1, and if
00061 *          JOBZ = 'V', LDZ >= max(1,N).
00062 *
00063 *  WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)
00064 *
00065 *  INFO    (output) INTEGER
00066 *          = 0:  successful exit.
00067 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00068 *          > 0:  if INFO = i, the algorithm failed to converge; i
00069 *                off-diagonal elements of an intermediate tridiagonal
00070 *                form did not converge to zero.
00071 *
00072 *  =====================================================================
00073 *
00074 *     .. Parameters ..
00075       DOUBLE PRECISION   ZERO, ONE
00076       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00077 *     ..
00078 *     .. Local Scalars ..
00079       LOGICAL            WANTZ
00080       INTEGER            IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE
00081       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
00082      $                   SMLNUM
00083 *     ..
00084 *     .. External Functions ..
00085       LOGICAL            LSAME
00086       DOUBLE PRECISION   DLAMCH, DLANSP
00087       EXTERNAL           LSAME, DLAMCH, DLANSP
00088 *     ..
00089 *     .. External Subroutines ..
00090       EXTERNAL           DOPGTR, DSCAL, DSPTRD, DSTEQR, DSTERF, XERBLA
00091 *     ..
00092 *     .. Intrinsic Functions ..
00093       INTRINSIC          SQRT
00094 *     ..
00095 *     .. Executable Statements ..
00096 *
00097 *     Test the input parameters.
00098 *
00099       WANTZ = LSAME( JOBZ, 'V' )
00100 *
00101       INFO = 0
00102       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
00103          INFO = -1
00104       ELSE IF( .NOT.( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) )
00105      $          THEN
00106          INFO = -2
00107       ELSE IF( N.LT.0 ) THEN
00108          INFO = -3
00109       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
00110          INFO = -7
00111       END IF
00112 *
00113       IF( INFO.NE.0 ) THEN
00114          CALL XERBLA( 'DSPEV ', -INFO )
00115          RETURN
00116       END IF
00117 *
00118 *     Quick return if possible
00119 *
00120       IF( N.EQ.0 )
00121      $   RETURN
00122 *
00123       IF( N.EQ.1 ) THEN
00124          W( 1 ) = AP( 1 )
00125          IF( WANTZ )
00126      $      Z( 1, 1 ) = ONE
00127          RETURN
00128       END IF
00129 *
00130 *     Get machine constants.
00131 *
00132       SAFMIN = DLAMCH( 'Safe minimum' )
00133       EPS = DLAMCH( 'Precision' )
00134       SMLNUM = SAFMIN / EPS
00135       BIGNUM = ONE / SMLNUM
00136       RMIN = SQRT( SMLNUM )
00137       RMAX = SQRT( BIGNUM )
00138 *
00139 *     Scale matrix to allowable range, if necessary.
00140 *
00141       ANRM = DLANSP( 'M', UPLO, N, AP, WORK )
00142       ISCALE = 0
00143       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
00144          ISCALE = 1
00145          SIGMA = RMIN / ANRM
00146       ELSE IF( ANRM.GT.RMAX ) THEN
00147          ISCALE = 1
00148          SIGMA = RMAX / ANRM
00149       END IF
00150       IF( ISCALE.EQ.1 ) THEN
00151          CALL DSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 )
00152       END IF
00153 *
00154 *     Call DSPTRD to reduce symmetric packed matrix to tridiagonal form.
00155 *
00156       INDE = 1
00157       INDTAU = INDE + N
00158       CALL DSPTRD( UPLO, N, AP, W, WORK( INDE ), WORK( INDTAU ), IINFO )
00159 *
00160 *     For eigenvalues only, call DSTERF.  For eigenvectors, first call
00161 *     DOPGTR to generate the orthogonal matrix, then call DSTEQR.
00162 *
00163       IF( .NOT.WANTZ ) THEN
00164          CALL DSTERF( N, W, WORK( INDE ), INFO )
00165       ELSE
00166          INDWRK = INDTAU + N
00167          CALL DOPGTR( UPLO, N, AP, WORK( INDTAU ), Z, LDZ,
00168      $                WORK( INDWRK ), IINFO )
00169          CALL DSTEQR( JOBZ, N, W, WORK( INDE ), Z, LDZ, WORK( INDTAU ),
00170      $                INFO )
00171       END IF
00172 *
00173 *     If matrix was scaled, then rescale eigenvalues appropriately.
00174 *
00175       IF( ISCALE.EQ.1 ) THEN
00176          IF( INFO.EQ.0 ) THEN
00177             IMAX = N
00178          ELSE
00179             IMAX = INFO - 1
00180          END IF
00181          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
00182       END IF
00183 *
00184       RETURN
00185 *
00186 *     End of DSPEV
00187 *
00188       END