LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sgeev.f
Go to the documentation of this file.
00001 *> \brief <b> SGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SGEEV + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgeev.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgeev.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgeev.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
00022 *                         LDVR, WORK, LWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          JOBVL, JOBVR
00026 *       INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       REAL               A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
00030 *      $                   WI( * ), WORK( * ), WR( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> SGEEV computes for an N-by-N real nonsymmetric matrix A, the
00040 *> eigenvalues and, optionally, the left and/or right eigenvectors.
00041 *>
00042 *> The right eigenvector v(j) of A satisfies
00043 *>                  A * v(j) = lambda(j) * v(j)
00044 *> where lambda(j) is its eigenvalue.
00045 *> The left eigenvector u(j) of A satisfies
00046 *>               u(j)**T * A = lambda(j) * u(j)**T
00047 *> where u(j)**T denotes the transpose of u(j).
00048 *>
00049 *> The computed eigenvectors are normalized to have Euclidean norm
00050 *> equal to 1 and largest component real.
00051 *> \endverbatim
00052 *
00053 *  Arguments:
00054 *  ==========
00055 *
00056 *> \param[in] JOBVL
00057 *> \verbatim
00058 *>          JOBVL is CHARACTER*1
00059 *>          = 'N': left eigenvectors of A are not computed;
00060 *>          = 'V': left eigenvectors of A are computed.
00061 *> \endverbatim
00062 *>
00063 *> \param[in] JOBVR
00064 *> \verbatim
00065 *>          JOBVR is CHARACTER*1
00066 *>          = 'N': right eigenvectors of A are not computed;
00067 *>          = 'V': right eigenvectors of A are computed.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] N
00071 *> \verbatim
00072 *>          N is INTEGER
00073 *>          The order of the matrix A. N >= 0.
00074 *> \endverbatim
00075 *>
00076 *> \param[in,out] A
00077 *> \verbatim
00078 *>          A is REAL array, dimension (LDA,N)
00079 *>          On entry, the N-by-N matrix A.
00080 *>          On exit, A has been overwritten.
00081 *> \endverbatim
00082 *>
00083 *> \param[in] LDA
00084 *> \verbatim
00085 *>          LDA is INTEGER
00086 *>          The leading dimension of the array A.  LDA >= max(1,N).
00087 *> \endverbatim
00088 *>
00089 *> \param[out] WR
00090 *> \verbatim
00091 *>          WR is REAL array, dimension (N)
00092 *> \endverbatim
00093 *>
00094 *> \param[out] WI
00095 *> \verbatim
00096 *>          WI is REAL array, dimension (N)
00097 *>          WR and WI contain the real and imaginary parts,
00098 *>          respectively, of the computed eigenvalues.  Complex
00099 *>          conjugate pairs of eigenvalues appear consecutively
00100 *>          with the eigenvalue having the positive imaginary part
00101 *>          first.
00102 *> \endverbatim
00103 *>
00104 *> \param[out] VL
00105 *> \verbatim
00106 *>          VL is REAL array, dimension (LDVL,N)
00107 *>          If JOBVL = 'V', the left eigenvectors u(j) are stored one
00108 *>          after another in the columns of VL, in the same order
00109 *>          as their eigenvalues.
00110 *>          If JOBVL = 'N', VL is not referenced.
00111 *>          If the j-th eigenvalue is real, then u(j) = VL(:,j),
00112 *>          the j-th column of VL.
00113 *>          If the j-th and (j+1)-st eigenvalues form a complex
00114 *>          conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
00115 *>          u(j+1) = VL(:,j) - i*VL(:,j+1).
00116 *> \endverbatim
00117 *>
00118 *> \param[in] LDVL
00119 *> \verbatim
00120 *>          LDVL is INTEGER
00121 *>          The leading dimension of the array VL.  LDVL >= 1; if
00122 *>          JOBVL = 'V', LDVL >= N.
00123 *> \endverbatim
00124 *>
00125 *> \param[out] VR
00126 *> \verbatim
00127 *>          VR is REAL array, dimension (LDVR,N)
00128 *>          If JOBVR = 'V', the right eigenvectors v(j) are stored one
00129 *>          after another in the columns of VR, in the same order
00130 *>          as their eigenvalues.
00131 *>          If JOBVR = 'N', VR is not referenced.
00132 *>          If the j-th eigenvalue is real, then v(j) = VR(:,j),
00133 *>          the j-th column of VR.
00134 *>          If the j-th and (j+1)-st eigenvalues form a complex
00135 *>          conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
00136 *>          v(j+1) = VR(:,j) - i*VR(:,j+1).
00137 *> \endverbatim
00138 *>
00139 *> \param[in] LDVR
00140 *> \verbatim
00141 *>          LDVR is INTEGER
00142 *>          The leading dimension of the array VR.  LDVR >= 1; if
00143 *>          JOBVR = 'V', LDVR >= N.
00144 *> \endverbatim
00145 *>
00146 *> \param[out] WORK
00147 *> \verbatim
00148 *>          WORK is REAL array, dimension (MAX(1,LWORK))
00149 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00150 *> \endverbatim
00151 *>
00152 *> \param[in] LWORK
00153 *> \verbatim
00154 *>          LWORK is INTEGER
00155 *>          The dimension of the array WORK.  LWORK >= max(1,3*N), and
00156 *>          if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N.  For good
00157 *>          performance, LWORK must generally be larger.
00158 *>
00159 *>          If LWORK = -1, then a workspace query is assumed; the routine
00160 *>          only calculates the optimal size of the WORK array, returns
00161 *>          this value as the first entry of the WORK array, and no error
00162 *>          message related to LWORK is issued by XERBLA.
00163 *> \endverbatim
00164 *>
00165 *> \param[out] INFO
00166 *> \verbatim
00167 *>          INFO is INTEGER
00168 *>          = 0:  successful exit
00169 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00170 *>          > 0:  if INFO = i, the QR algorithm failed to compute all the
00171 *>                eigenvalues, and no eigenvectors have been computed;
00172 *>                elements i+1:N of WR and WI contain eigenvalues which
00173 *>                have converged.
00174 *> \endverbatim
00175 *
00176 *  Authors:
00177 *  ========
00178 *
00179 *> \author Univ. of Tennessee 
00180 *> \author Univ. of California Berkeley 
00181 *> \author Univ. of Colorado Denver 
00182 *> \author NAG Ltd. 
00183 *
00184 *> \date November 2011
00185 *
00186 *> \ingroup realGEeigen
00187 *
00188 *  =====================================================================
00189       SUBROUTINE SGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
00190      $                  LDVR, WORK, LWORK, INFO )
00191 *
00192 *  -- LAPACK driver routine (version 3.4.0) --
00193 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00194 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00195 *     November 2011
00196 *
00197 *     .. Scalar Arguments ..
00198       CHARACTER          JOBVL, JOBVR
00199       INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
00200 *     ..
00201 *     .. Array Arguments ..
00202       REAL               A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
00203      $                   WI( * ), WORK( * ), WR( * )
00204 *     ..
00205 *
00206 *  =====================================================================
00207 *
00208 *     .. Parameters ..
00209       REAL               ZERO, ONE
00210       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00211 *     ..
00212 *     .. Local Scalars ..
00213       LOGICAL            LQUERY, SCALEA, WANTVL, WANTVR
00214       CHARACTER          SIDE
00215       INTEGER            HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
00216      $                   MAXWRK, MINWRK, NOUT
00217       REAL               ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
00218      $                   SN
00219 *     ..
00220 *     .. Local Arrays ..
00221       LOGICAL            SELECT( 1 )
00222       REAL               DUM( 1 )
00223 *     ..
00224 *     .. External Subroutines ..
00225       EXTERNAL           SGEBAK, SGEBAL, SGEHRD, SHSEQR, SLABAD, SLACPY,
00226      $                   SLARTG, SLASCL, SORGHR, SROT, SSCAL, STREVC,
00227      $                   XERBLA
00228 *     ..
00229 *     .. External Functions ..
00230       LOGICAL            LSAME
00231       INTEGER            ILAENV, ISAMAX
00232       REAL               SLAMCH, SLANGE, SLAPY2, SNRM2
00233       EXTERNAL           LSAME, ILAENV, ISAMAX, SLAMCH, SLANGE, SLAPY2,
00234      $                   SNRM2
00235 *     ..
00236 *     .. Intrinsic Functions ..
00237       INTRINSIC          MAX, SQRT
00238 *     ..
00239 *     .. Executable Statements ..
00240 *
00241 *     Test the input arguments
00242 *
00243       INFO = 0
00244       LQUERY = ( LWORK.EQ.-1 )
00245       WANTVL = LSAME( JOBVL, 'V' )
00246       WANTVR = LSAME( JOBVR, 'V' )
00247       IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
00248          INFO = -1
00249       ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
00250          INFO = -2
00251       ELSE IF( N.LT.0 ) THEN
00252          INFO = -3
00253       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00254          INFO = -5
00255       ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
00256          INFO = -9
00257       ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
00258          INFO = -11
00259       END IF
00260 *
00261 *     Compute workspace
00262 *      (Note: Comments in the code beginning "Workspace:" describe the
00263 *       minimal amount of workspace needed at that point in the code,
00264 *       as well as the preferred amount for good performance.
00265 *       NB refers to the optimal block size for the immediately
00266 *       following subroutine, as returned by ILAENV.
00267 *       HSWORK refers to the workspace preferred by SHSEQR, as
00268 *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
00269 *       the worst case.)
00270 *
00271       IF( INFO.EQ.0 ) THEN
00272          IF( N.EQ.0 ) THEN
00273             MINWRK = 1
00274             MAXWRK = 1
00275          ELSE
00276             MAXWRK = 2*N + N*ILAENV( 1, 'SGEHRD', ' ', N, 1, N, 0 )
00277             IF( WANTVL ) THEN
00278                MINWRK = 4*N
00279                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
00280      $                       'SORGHR', ' ', N, 1, N, -1 ) )
00281                CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
00282      $                WORK, -1, INFO )
00283                HSWORK = WORK( 1 )
00284                MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
00285                MAXWRK = MAX( MAXWRK, 4*N )
00286             ELSE IF( WANTVR ) THEN
00287                MINWRK = 4*N
00288                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
00289      $                       'SORGHR', ' ', N, 1, N, -1 ) )
00290                CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
00291      $                WORK, -1, INFO )
00292                HSWORK = WORK( 1 )
00293                MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
00294                MAXWRK = MAX( MAXWRK, 4*N )
00295             ELSE 
00296                MINWRK = 3*N
00297                CALL SHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR,
00298      $                WORK, -1, INFO )
00299                HSWORK = WORK( 1 )
00300                MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
00301             END IF
00302             MAXWRK = MAX( MAXWRK, MINWRK )
00303          END IF
00304          WORK( 1 ) = MAXWRK
00305 *
00306          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
00307             INFO = -13
00308          END IF
00309       END IF
00310 *
00311       IF( INFO.NE.0 ) THEN
00312          CALL XERBLA( 'SGEEV ', -INFO )
00313          RETURN
00314       ELSE IF( LQUERY ) THEN
00315          RETURN
00316       END IF
00317 *
00318 *     Quick return if possible
00319 *
00320       IF( N.EQ.0 )
00321      $   RETURN
00322 *
00323 *     Get machine constants
00324 *
00325       EPS = SLAMCH( 'P' )
00326       SMLNUM = SLAMCH( 'S' )
00327       BIGNUM = ONE / SMLNUM
00328       CALL SLABAD( SMLNUM, BIGNUM )
00329       SMLNUM = SQRT( SMLNUM ) / EPS
00330       BIGNUM = ONE / SMLNUM
00331 *
00332 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00333 *
00334       ANRM = SLANGE( 'M', N, N, A, LDA, DUM )
00335       SCALEA = .FALSE.
00336       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00337          SCALEA = .TRUE.
00338          CSCALE = SMLNUM
00339       ELSE IF( ANRM.GT.BIGNUM ) THEN
00340          SCALEA = .TRUE.
00341          CSCALE = BIGNUM
00342       END IF
00343       IF( SCALEA )
00344      $   CALL SLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
00345 *
00346 *     Balance the matrix
00347 *     (Workspace: need N)
00348 *
00349       IBAL = 1
00350       CALL SGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
00351 *
00352 *     Reduce to upper Hessenberg form
00353 *     (Workspace: need 3*N, prefer 2*N+N*NB)
00354 *
00355       ITAU = IBAL + N
00356       IWRK = ITAU + N
00357       CALL SGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
00358      $             LWORK-IWRK+1, IERR )
00359 *
00360       IF( WANTVL ) THEN
00361 *
00362 *        Want left eigenvectors
00363 *        Copy Householder vectors to VL
00364 *
00365          SIDE = 'L'
00366          CALL SLACPY( 'L', N, N, A, LDA, VL, LDVL )
00367 *
00368 *        Generate orthogonal matrix in VL
00369 *        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
00370 *
00371          CALL SORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
00372      $                LWORK-IWRK+1, IERR )
00373 *
00374 *        Perform QR iteration, accumulating Schur vectors in VL
00375 *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
00376 *
00377          IWRK = ITAU
00378          CALL SHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
00379      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00380 *
00381          IF( WANTVR ) THEN
00382 *
00383 *           Want left and right eigenvectors
00384 *           Copy Schur vectors to VR
00385 *
00386             SIDE = 'B'
00387             CALL SLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
00388          END IF
00389 *
00390       ELSE IF( WANTVR ) THEN
00391 *
00392 *        Want right eigenvectors
00393 *        Copy Householder vectors to VR
00394 *
00395          SIDE = 'R'
00396          CALL SLACPY( 'L', N, N, A, LDA, VR, LDVR )
00397 *
00398 *        Generate orthogonal matrix in VR
00399 *        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
00400 *
00401          CALL SORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
00402      $                LWORK-IWRK+1, IERR )
00403 *
00404 *        Perform QR iteration, accumulating Schur vectors in VR
00405 *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
00406 *
00407          IWRK = ITAU
00408          CALL SHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
00409      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00410 *
00411       ELSE
00412 *
00413 *        Compute eigenvalues only
00414 *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
00415 *
00416          IWRK = ITAU
00417          CALL SHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
00418      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00419       END IF
00420 *
00421 *     If INFO > 0 from SHSEQR, then quit
00422 *
00423       IF( INFO.GT.0 )
00424      $   GO TO 50
00425 *
00426       IF( WANTVL .OR. WANTVR ) THEN
00427 *
00428 *        Compute left and/or right eigenvectors
00429 *        (Workspace: need 4*N)
00430 *
00431          CALL STREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
00432      $                N, NOUT, WORK( IWRK ), IERR )
00433       END IF
00434 *
00435       IF( WANTVL ) THEN
00436 *
00437 *        Undo balancing of left eigenvectors
00438 *        (Workspace: need N)
00439 *
00440          CALL SGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,
00441      $                IERR )
00442 *
00443 *        Normalize left eigenvectors and make largest component real
00444 *
00445          DO 20 I = 1, N
00446             IF( WI( I ).EQ.ZERO ) THEN
00447                SCL = ONE / SNRM2( N, VL( 1, I ), 1 )
00448                CALL SSCAL( N, SCL, VL( 1, I ), 1 )
00449             ELSE IF( WI( I ).GT.ZERO ) THEN
00450                SCL = ONE / SLAPY2( SNRM2( N, VL( 1, I ), 1 ),
00451      $               SNRM2( N, VL( 1, I+1 ), 1 ) )
00452                CALL SSCAL( N, SCL, VL( 1, I ), 1 )
00453                CALL SSCAL( N, SCL, VL( 1, I+1 ), 1 )
00454                DO 10 K = 1, N
00455                   WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2
00456    10          CONTINUE
00457                K = ISAMAX( N, WORK( IWRK ), 1 )
00458                CALL SLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
00459                CALL SROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
00460                VL( K, I+1 ) = ZERO
00461             END IF
00462    20    CONTINUE
00463       END IF
00464 *
00465       IF( WANTVR ) THEN
00466 *
00467 *        Undo balancing of right eigenvectors
00468 *        (Workspace: need N)
00469 *
00470          CALL SGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,
00471      $                IERR )
00472 *
00473 *        Normalize right eigenvectors and make largest component real
00474 *
00475          DO 40 I = 1, N
00476             IF( WI( I ).EQ.ZERO ) THEN
00477                SCL = ONE / SNRM2( N, VR( 1, I ), 1 )
00478                CALL SSCAL( N, SCL, VR( 1, I ), 1 )
00479             ELSE IF( WI( I ).GT.ZERO ) THEN
00480                SCL = ONE / SLAPY2( SNRM2( N, VR( 1, I ), 1 ),
00481      $               SNRM2( N, VR( 1, I+1 ), 1 ) )
00482                CALL SSCAL( N, SCL, VR( 1, I ), 1 )
00483                CALL SSCAL( N, SCL, VR( 1, I+1 ), 1 )
00484                DO 30 K = 1, N
00485                   WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2
00486    30          CONTINUE
00487                K = ISAMAX( N, WORK( IWRK ), 1 )
00488                CALL SLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
00489                CALL SROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
00490                VR( K, I+1 ) = ZERO
00491             END IF
00492    40    CONTINUE
00493       END IF
00494 *
00495 *     Undo scaling if necessary
00496 *
00497    50 CONTINUE
00498       IF( SCALEA ) THEN
00499          CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
00500      $                MAX( N-INFO, 1 ), IERR )
00501          CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
00502      $                MAX( N-INFO, 1 ), IERR )
00503          IF( INFO.GT.0 ) THEN
00504             CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
00505      $                   IERR )
00506             CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
00507      $                   IERR )
00508          END IF
00509       END IF
00510 *
00511       WORK( 1 ) = MAXWRK
00512       RETURN
00513 *
00514 *     End of SGEEV
00515 *
00516       END
 All Files Functions