LAPACK: SRC/sgees.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, WR, WI,
00002      $                  VS, LDVS, WORK, LWORK, BWORK, 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          JOBVS, SORT
00011       INTEGER            INFO, LDA, LDVS, LWORK, N, SDIM
00012 *     ..
00013 *     .. Array Arguments ..
00014       LOGICAL            BWORK( * )
00015       REAL               A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
00016      $                   WR( * )
00017 *     ..
00018 *     .. Function Arguments ..
00019       LOGICAL            SELECT
00020       EXTERNAL           SELECT
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  SGEES computes for an N-by-N real nonsymmetric matrix A, the
00027 *  eigenvalues, the real Schur form T, and, optionally, the matrix of
00028 *  Schur vectors Z.  This gives the Schur factorization A = Z*T*(Z**T).
00029 *
00030 *  Optionally, it also orders the eigenvalues on the diagonal of the
00031 *  real Schur form so that selected eigenvalues are at the top left.
00032 *  The leading columns of Z then form an orthonormal basis for the
00033 *  invariant subspace corresponding to the selected eigenvalues.
00034 *
00035 *  A matrix is in real Schur form if it is upper quasi-triangular with
00036 *  1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in the
00037 *  form
00038 *          [  a  b  ]
00039 *          [  c  a  ]
00040 *
00041 *  where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
00042 *
00043 *  Arguments
00044 *  =========
00045 *
00046 *  JOBVS   (input) CHARACTER*1
00047 *          = 'N': Schur vectors are not computed;
00048 *          = 'V': Schur vectors are computed.
00049 *
00050 *  SORT    (input) CHARACTER*1
00051 *          Specifies whether or not to order the eigenvalues on the
00052 *          diagonal of the Schur form.
00053 *          = 'N': Eigenvalues are not ordered;
00054 *          = 'S': Eigenvalues are ordered (see SELECT).
00055 *
00056 *  SELECT  (external procedure) LOGICAL FUNCTION of two REAL arguments
00057 *          SELECT must be declared EXTERNAL in the calling subroutine.
00058 *          If SORT = 'S', SELECT is used to select eigenvalues to sort
00059 *          to the top left of the Schur form.
00060 *          If SORT = 'N', SELECT is not referenced.
00061 *          An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
00062 *          SELECT(WR(j),WI(j)) is true; i.e., if either one of a complex
00063 *          conjugate pair of eigenvalues is selected, then both complex
00064 *          eigenvalues are selected.
00065 *          Note that a selected complex eigenvalue may no longer
00066 *          satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
00067 *          ordering may change the value of complex eigenvalues
00068 *          (especially if the eigenvalue is ill-conditioned); in this
00069 *          case INFO is set to N+2 (see INFO below).
00070 *
00071 *  N       (input) INTEGER
00072 *          The order of the matrix A. N >= 0.
00073 *
00074 *  A       (input/output) REAL array, dimension (LDA,N)
00075 *          On entry, the N-by-N matrix A.
00076 *          On exit, A has been overwritten by its real Schur form T.
00077 *
00078 *  LDA     (input) INTEGER
00079 *          The leading dimension of the array A.  LDA >= max(1,N).
00080 *
00081 *  SDIM    (output) INTEGER
00082 *          If SORT = 'N', SDIM = 0.
00083 *          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
00084 *                         for which SELECT is true. (Complex conjugate
00085 *                         pairs for which SELECT is true for either
00086 *                         eigenvalue count as 2.)
00087 *
00088 *  WR      (output) REAL array, dimension (N)
00089 *  WI      (output) REAL array, dimension (N)
00090 *          WR and WI contain the real and imaginary parts,
00091 *          respectively, of the computed eigenvalues in the same order
00092 *          that they appear on the diagonal of the output Schur form T.
00093 *          Complex conjugate pairs of eigenvalues will appear
00094 *          consecutively with the eigenvalue having the positive
00095 *          imaginary part first.
00096 *
00097 *  VS      (output) REAL array, dimension (LDVS,N)
00098 *          If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
00099 *          vectors.
00100 *          If JOBVS = 'N', VS is not referenced.
00101 *
00102 *  LDVS    (input) INTEGER
00103 *          The leading dimension of the array VS.  LDVS >= 1; if
00104 *          JOBVS = 'V', LDVS >= N.
00105 *
00106 *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
00107 *          On exit, if INFO = 0, WORK(1) contains the optimal LWORK.
00108 *
00109 *  LWORK   (input) INTEGER
00110 *          The dimension of the array WORK.  LWORK >= max(1,3*N).
00111 *          For good performance, LWORK must generally be larger.
00112 *
00113 *          If LWORK = -1, then a workspace query is assumed; the routine
00114 *          only calculates the optimal size of the WORK array, returns
00115 *          this value as the first entry of the WORK array, and no error
00116 *          message related to LWORK is issued by XERBLA.
00117 *
00118 *  BWORK   (workspace) LOGICAL array, dimension (N)
00119 *          Not referenced if SORT = 'N'.
00120 *
00121 *  INFO    (output) INTEGER
00122 *          = 0: successful exit
00123 *          < 0: if INFO = -i, the i-th argument had an illegal value.
00124 *          > 0: if INFO = i, and i is
00125 *             <= N: the QR algorithm failed to compute all the
00126 *                   eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
00127 *                   contain those eigenvalues which have converged; if
00128 *                   JOBVS = 'V', VS contains the matrix which reduces A
00129 *                   to its partially converged Schur form.
00130 *             = N+1: the eigenvalues could not be reordered because some
00131 *                   eigenvalues were too close to separate (the problem
00132 *                   is very ill-conditioned);
00133 *             = N+2: after reordering, roundoff changed values of some
00134 *                   complex eigenvalues so that leading eigenvalues in
00135 *                   the Schur form no longer satisfy SELECT=.TRUE.  This
00136 *                   could also be caused by underflow due to scaling.
00137 *
00138 *  =====================================================================
00139 *
00140 *     .. Parameters ..
00141       REAL               ZERO, ONE
00142       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00143 *     ..
00144 *     .. Local Scalars ..
00145       LOGICAL            CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTST,
00146      $                   WANTVS
00147       INTEGER            HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
00148      $                   IHI, ILO, INXT, IP, ITAU, IWRK, MAXWRK, MINWRK
00149       REAL               ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
00150 *     ..
00151 *     .. Local Arrays ..
00152       INTEGER            IDUM( 1 )
00153       REAL               DUM( 1 )
00154 *     ..
00155 *     .. External Subroutines ..
00156       EXTERNAL           SCOPY, SGEBAK, SGEBAL, SGEHRD, SHSEQR, SLABAD,
00157      $                   SLACPY, SLASCL, SORGHR, SSWAP, STRSEN, XERBLA
00158 *     ..
00159 *     .. External Functions ..
00160       LOGICAL            LSAME
00161       INTEGER            ILAENV
00162       REAL               SLAMCH, SLANGE
00163       EXTERNAL           LSAME, ILAENV, SLAMCH, SLANGE
00164 *     ..
00165 *     .. Intrinsic Functions ..
00166       INTRINSIC          MAX, SQRT
00167 *     ..
00168 *     .. Executable Statements ..
00169 *
00170 *     Test the input arguments
00171 *
00172       INFO = 0
00173       LQUERY = ( LWORK.EQ.-1 )
00174       WANTVS = LSAME( JOBVS, 'V' )
00175       WANTST = LSAME( SORT, 'S' )
00176       IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
00177          INFO = -1
00178       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
00179          INFO = -2
00180       ELSE IF( N.LT.0 ) THEN
00181          INFO = -4
00182       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00183          INFO = -6
00184       ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
00185          INFO = -11
00186       END IF
00187 *
00188 *     Compute workspace
00189 *      (Note: Comments in the code beginning "Workspace:" describe the
00190 *       minimal amount of workspace needed at that point in the code,
00191 *       as well as the preferred amount for good performance.
00192 *       NB refers to the optimal block size for the immediately
00193 *       following subroutine, as returned by ILAENV.
00194 *       HSWORK refers to the workspace preferred by SHSEQR, as
00195 *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
00196 *       the worst case.)
00197 *
00198       IF( INFO.EQ.0 ) THEN
00199          IF( N.EQ.0 ) THEN
00200             MINWRK = 1
00201             MAXWRK = 1
00202          ELSE
00203             MAXWRK = 2*N + N*ILAENV( 1, 'SGEHRD', ' ', N, 1, N, 0 )
00204             MINWRK = 3*N
00205 *
00206             CALL SHSEQR( 'S', JOBVS, N, 1, N, A, LDA, WR, WI, VS, LDVS,
00207      $             WORK, -1, IEVAL )
00208             HSWORK = WORK( 1 )
00209 *
00210             IF( .NOT.WANTVS ) THEN
00211                MAXWRK = MAX( MAXWRK, N + HSWORK )
00212             ELSE
00213                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
00214      $                       'SORGHR', ' ', N, 1, N, -1 ) )
00215                MAXWRK = MAX( MAXWRK, N + HSWORK )
00216             END IF
00217          END IF
00218          WORK( 1 ) = MAXWRK
00219 *
00220          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
00221             INFO = -13
00222          END IF
00223       END IF
00224 *
00225       IF( INFO.NE.0 ) THEN
00226          CALL XERBLA( 'SGEES ', -INFO )
00227          RETURN
00228       ELSE IF( LQUERY ) THEN
00229          RETURN
00230       END IF
00231 *
00232 *     Quick return if possible
00233 *
00234       IF( N.EQ.0 ) THEN
00235          SDIM = 0
00236          RETURN
00237       END IF
00238 *
00239 *     Get machine constants
00240 *
00241       EPS = SLAMCH( 'P' )
00242       SMLNUM = SLAMCH( 'S' )
00243       BIGNUM = ONE / SMLNUM
00244       CALL SLABAD( SMLNUM, BIGNUM )
00245       SMLNUM = SQRT( SMLNUM ) / EPS
00246       BIGNUM = ONE / SMLNUM
00247 *
00248 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00249 *
00250       ANRM = SLANGE( 'M', N, N, A, LDA, DUM )
00251       SCALEA = .FALSE.
00252       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00253          SCALEA = .TRUE.
00254          CSCALE = SMLNUM
00255       ELSE IF( ANRM.GT.BIGNUM ) THEN
00256          SCALEA = .TRUE.
00257          CSCALE = BIGNUM
00258       END IF
00259       IF( SCALEA )
00260      $   CALL SLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
00261 *
00262 *     Permute the matrix to make it more nearly triangular
00263 *     (Workspace: need N)
00264 *
00265       IBAL = 1
00266       CALL SGEBAL( 'P', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
00267 *
00268 *     Reduce to upper Hessenberg form
00269 *     (Workspace: need 3*N, prefer 2*N+N*NB)
00270 *
00271       ITAU = N + IBAL
00272       IWRK = N + ITAU
00273       CALL SGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
00274      $             LWORK-IWRK+1, IERR )
00275 *
00276       IF( WANTVS ) THEN
00277 *
00278 *        Copy Householder vectors to VS
00279 *
00280          CALL SLACPY( 'L', N, N, A, LDA, VS, LDVS )
00281 *
00282 *        Generate orthogonal matrix in VS
00283 *        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
00284 *
00285          CALL SORGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
00286      $                LWORK-IWRK+1, IERR )
00287       END IF
00288 *
00289       SDIM = 0
00290 *
00291 *     Perform QR iteration, accumulating Schur vectors in VS if desired
00292 *     (Workspace: need N+1, prefer N+HSWORK (see comments) )
00293 *
00294       IWRK = ITAU
00295       CALL SHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, WR, WI, VS, LDVS,
00296      $             WORK( IWRK ), LWORK-IWRK+1, IEVAL )
00297       IF( IEVAL.GT.0 )
00298      $   INFO = IEVAL
00299 *
00300 *     Sort eigenvalues if desired
00301 *
00302       IF( WANTST .AND. INFO.EQ.0 ) THEN
00303          IF( SCALEA ) THEN
00304             CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WR, N, IERR )
00305             CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WI, N, IERR )
00306          END IF
00307          DO 10 I = 1, N
00308             BWORK( I ) = SELECT( WR( I ), WI( I ) )
00309    10    CONTINUE
00310 *
00311 *        Reorder eigenvalues and transform Schur vectors
00312 *        (Workspace: none needed)
00313 *
00314          CALL STRSEN( 'N', JOBVS, BWORK, N, A, LDA, VS, LDVS, WR, WI,
00315      $                SDIM, S, SEP, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1,
00316      $                ICOND )
00317          IF( ICOND.GT.0 )
00318      $      INFO = N + ICOND
00319       END IF
00320 *
00321       IF( WANTVS ) THEN
00322 *
00323 *        Undo balancing
00324 *        (Workspace: need N)
00325 *
00326          CALL SGEBAK( 'P', 'R', N, ILO, IHI, WORK( IBAL ), N, VS, LDVS,
00327      $                IERR )
00328       END IF
00329 *
00330       IF( SCALEA ) THEN
00331 *
00332 *        Undo scaling for the Schur form of A
00333 *
00334          CALL SLASCL( 'H', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
00335          CALL SCOPY( N, A, LDA+1, WR, 1 )
00336          IF( CSCALE.EQ.SMLNUM ) THEN
00337 *
00338 *           If scaling back towards underflow, adjust WI if an
00339 *           offdiagonal element of a 2-by-2 block in the Schur form
00340 *           underflows.
00341 *
00342             IF( IEVAL.GT.0 ) THEN
00343                I1 = IEVAL + 1
00344                I2 = IHI - 1
00345                CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI,
00346      $                      MAX( ILO-1, 1 ), IERR )
00347             ELSE IF( WANTST ) THEN
00348                I1 = 1
00349                I2 = N - 1
00350             ELSE
00351                I1 = ILO
00352                I2 = IHI - 1
00353             END IF
00354             INXT = I1 - 1
00355             DO 20 I = I1, I2
00356                IF( I.LT.INXT )
00357      $            GO TO 20
00358                IF( WI( I ).EQ.ZERO ) THEN
00359                   INXT = I + 1
00360                ELSE
00361                   IF( A( I+1, I ).EQ.ZERO ) THEN
00362                      WI( I ) = ZERO
00363                      WI( I+1 ) = ZERO
00364                   ELSE IF( A( I+1, I ).NE.ZERO .AND. A( I, I+1 ).EQ.
00365      $                     ZERO ) THEN
00366                      WI( I ) = ZERO
00367                      WI( I+1 ) = ZERO
00368                      IF( I.GT.1 )
00369      $                  CALL SSWAP( I-1, A( 1, I ), 1, A( 1, I+1 ), 1 )
00370                      IF( N.GT.I+1 )
00371      $                  CALL SSWAP( N-I-1, A( I, I+2 ), LDA,
00372      $                              A( I+1, I+2 ), LDA )
00373                      IF( WANTVS ) THEN
00374                         CALL SSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
00375                      END IF
00376                      A( I, I+1 ) = A( I+1, I )
00377                      A( I+1, I ) = ZERO
00378                   END IF
00379                   INXT = I + 2
00380                END IF
00381    20       CONTINUE
00382          END IF
00383 *
00384 *        Undo scaling for the imaginary part of the eigenvalues
00385 *
00386          CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N-IEVAL, 1,
00387      $                WI( IEVAL+1 ), MAX( N-IEVAL, 1 ), IERR )
00388       END IF
00389 *
00390       IF( WANTST .AND. INFO.EQ.0 ) THEN
00391 *
00392 *        Check if reordering successful
00393 *
00394          LASTSL = .TRUE.
00395          LST2SL = .TRUE.
00396          SDIM = 0
00397          IP = 0
00398          DO 30 I = 1, N
00399             CURSL = SELECT( WR( I ), WI( I ) )
00400             IF( WI( I ).EQ.ZERO ) THEN
00401                IF( CURSL )
00402      $            SDIM = SDIM + 1
00403                IP = 0
00404                IF( CURSL .AND. .NOT.LASTSL )
00405      $            INFO = N + 2
00406             ELSE
00407                IF( IP.EQ.1 ) THEN
00408 *
00409 *                 Last eigenvalue of conjugate pair
00410 *
00411                   CURSL = CURSL .OR. LASTSL
00412                   LASTSL = CURSL
00413                   IF( CURSL )
00414      $               SDIM = SDIM + 2
00415                   IP = -1
00416                   IF( CURSL .AND. .NOT.LST2SL )
00417      $               INFO = N + 2
00418                ELSE
00419 *
00420 *                 First eigenvalue of conjugate pair
00421 *
00422                   IP = 1
00423                END IF
00424             END IF
00425             LST2SL = LASTSL
00426             LASTSL = CURSL
00427    30    CONTINUE
00428       END IF
00429 *
00430       WORK( 1 ) = MAXWRK
00431       RETURN
00432 *
00433 *     End of SGEES
00434 *
00435       END