LAPACK 3.3.0

cgees.f

Go to the documentation of this file.
00001       SUBROUTINE CGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
00002      $                  LDVS, WORK, LWORK, RWORK, 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               RWORK( * )
00016       COMPLEX            A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
00017 *     ..
00018 *     .. Function Arguments ..
00019       LOGICAL            SELECT
00020       EXTERNAL           SELECT
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  CGEES computes for an N-by-N complex nonsymmetric matrix A, the
00027 *  eigenvalues, the Schur form T, and, optionally, the matrix of Schur
00028 *  vectors Z.  This gives the Schur factorization A = Z*T*(Z**H).
00029 *
00030 *  Optionally, it also orders the eigenvalues on the diagonal of the
00031 *  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 complex matrix is in Schur form if it is upper triangular.
00036 *
00037 *  Arguments
00038 *  =========
00039 *
00040 *  JOBVS   (input) CHARACTER*1
00041 *          = 'N': Schur vectors are not computed;
00042 *          = 'V': Schur vectors are computed.
00043 *
00044 *  SORT    (input) CHARACTER*1
00045 *          Specifies whether or not to order the eigenvalues on the
00046 *          diagonal of the Schur form.
00047 *          = 'N': Eigenvalues are not ordered:
00048 *          = 'S': Eigenvalues are ordered (see SELECT).
00049 *
00050 *  SELECT  (external procedure) LOGICAL FUNCTION of one COMPLEX argument
00051 *          SELECT must be declared EXTERNAL in the calling subroutine.
00052 *          If SORT = 'S', SELECT is used to select eigenvalues to order
00053 *          to the top left of the Schur form.
00054 *          IF SORT = 'N', SELECT is not referenced.
00055 *          The eigenvalue W(j) is selected if SELECT(W(j)) is true.
00056 *
00057 *  N       (input) INTEGER
00058 *          The order of the matrix A. N >= 0.
00059 *
00060 *  A       (input/output) COMPLEX array, dimension (LDA,N)
00061 *          On entry, the N-by-N matrix A.
00062 *          On exit, A has been overwritten by its Schur form T.
00063 *
00064 *  LDA     (input) INTEGER
00065 *          The leading dimension of the array A.  LDA >= max(1,N).
00066 *
00067 *  SDIM    (output) INTEGER
00068 *          If SORT = 'N', SDIM = 0.
00069 *          If SORT = 'S', SDIM = number of eigenvalues for which
00070 *                         SELECT is true.
00071 *
00072 *  W       (output) COMPLEX array, dimension (N)
00073 *          W contains the computed eigenvalues, in the same order that
00074 *          they appear on the diagonal of the output Schur form T.
00075 *
00076 *  VS      (output) COMPLEX array, dimension (LDVS,N)
00077 *          If JOBVS = 'V', VS contains the unitary matrix Z of Schur
00078 *          vectors.
00079 *          If JOBVS = 'N', VS is not referenced.
00080 *
00081 *  LDVS    (input) INTEGER
00082 *          The leading dimension of the array VS.  LDVS >= 1; if
00083 *          JOBVS = 'V', LDVS >= N.
00084 *
00085 *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
00086 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00087 *
00088 *  LWORK   (input) INTEGER
00089 *          The dimension of the array WORK.  LWORK >= max(1,2*N).
00090 *          For good performance, LWORK must generally be larger.
00091 *
00092 *          If LWORK = -1, then a workspace query is assumed; the routine
00093 *          only calculates the optimal size of the WORK array, returns
00094 *          this value as the first entry of the WORK array, and no error
00095 *          message related to LWORK is issued by XERBLA.
00096 *
00097 *  RWORK   (workspace) REAL array, dimension (N)
00098 *
00099 *  BWORK   (workspace) LOGICAL array, dimension (N)
00100 *          Not referenced if SORT = 'N'.
00101 *
00102 *  INFO    (output) INTEGER
00103 *          = 0: successful exit
00104 *          < 0: if INFO = -i, the i-th argument had an illegal value.
00105 *          > 0: if INFO = i, and i is
00106 *               <= N:  the QR algorithm failed to compute all the
00107 *                      eigenvalues; elements 1:ILO-1 and i+1:N of W
00108 *                      contain those eigenvalues which have converged;
00109 *                      if JOBVS = 'V', VS contains the matrix which
00110 *                      reduces A to its partially converged Schur form.
00111 *               = N+1: the eigenvalues could not be reordered because
00112 *                      some eigenvalues were too close to separate (the
00113 *                      problem is very ill-conditioned);
00114 *               = N+2: after reordering, roundoff changed values of
00115 *                      some complex eigenvalues so that leading
00116 *                      eigenvalues in the Schur form no longer satisfy
00117 *                      SELECT = .TRUE..  This could also be caused by
00118 *                      underflow due to scaling.
00119 *
00120 *  =====================================================================
00121 *
00122 *     .. Parameters ..
00123       REAL               ZERO, ONE
00124       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00125 *     ..
00126 *     .. Local Scalars ..
00127       LOGICAL            LQUERY, SCALEA, WANTST, WANTVS
00128       INTEGER            HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
00129      $                   ITAU, IWRK, MAXWRK, MINWRK
00130       REAL               ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
00131 *     ..
00132 *     .. Local Arrays ..
00133       REAL               DUM( 1 )
00134 *     ..
00135 *     .. External Subroutines ..
00136       EXTERNAL           CCOPY, CGEBAK, CGEBAL, CGEHRD, CHSEQR, CLACPY,
00137      $                   CLASCL, CTRSEN, CUNGHR, SLABAD, XERBLA
00138 *     ..
00139 *     .. External Functions ..
00140       LOGICAL            LSAME
00141       INTEGER            ILAENV
00142       REAL               CLANGE, SLAMCH
00143       EXTERNAL           LSAME, ILAENV, CLANGE, SLAMCH
00144 *     ..
00145 *     .. Intrinsic Functions ..
00146       INTRINSIC          MAX, SQRT
00147 *     ..
00148 *     .. Executable Statements ..
00149 *
00150 *     Test the input arguments
00151 *
00152       INFO = 0
00153       LQUERY = ( LWORK.EQ.-1 )
00154       WANTVS = LSAME( JOBVS, 'V' )
00155       WANTST = LSAME( SORT, 'S' )
00156       IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
00157          INFO = -1
00158       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
00159          INFO = -2
00160       ELSE IF( N.LT.0 ) THEN
00161          INFO = -4
00162       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00163          INFO = -6
00164       ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
00165          INFO = -10
00166       END IF
00167 *
00168 *     Compute workspace
00169 *      (Note: Comments in the code beginning "Workspace:" describe the
00170 *       minimal amount of workspace needed at that point in the code,
00171 *       as well as the preferred amount for good performance.
00172 *       CWorkspace refers to complex workspace, and RWorkspace to real
00173 *       workspace. NB refers to the optimal block size for the
00174 *       immediately following subroutine, as returned by ILAENV.
00175 *       HSWORK refers to the workspace preferred by CHSEQR, as
00176 *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
00177 *       the worst case.)
00178 *
00179       IF( INFO.EQ.0 ) THEN
00180          IF( N.EQ.0 ) THEN
00181             MINWRK = 1
00182             MAXWRK = 1
00183          ELSE
00184             MAXWRK = N + N*ILAENV( 1, 'CGEHRD', ' ', N, 1, N, 0 )
00185             MINWRK = 2*N
00186 *
00187             CALL CHSEQR( 'S', JOBVS, N, 1, N, A, LDA, W, VS, LDVS,
00188      $             WORK, -1, IEVAL )
00189             HSWORK = WORK( 1 )
00190 *
00191             IF( .NOT.WANTVS ) THEN
00192                MAXWRK = MAX( MAXWRK, HSWORK )
00193             ELSE
00194                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'CUNGHR',
00195      $                       ' ', N, 1, N, -1 ) )
00196                MAXWRK = MAX( MAXWRK, HSWORK )
00197             END IF
00198          END IF
00199          WORK( 1 ) = MAXWRK
00200 *
00201          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
00202             INFO = -12
00203          END IF
00204       END IF
00205 *
00206       IF( INFO.NE.0 ) THEN
00207          CALL XERBLA( 'CGEES ', -INFO )
00208          RETURN
00209       ELSE IF( LQUERY ) THEN
00210          RETURN
00211       END IF
00212 *
00213 *     Quick return if possible
00214 *
00215       IF( N.EQ.0 ) THEN
00216          SDIM = 0
00217          RETURN
00218       END IF
00219 *
00220 *     Get machine constants
00221 *
00222       EPS = SLAMCH( 'P' )
00223       SMLNUM = SLAMCH( 'S' )
00224       BIGNUM = ONE / SMLNUM
00225       CALL SLABAD( SMLNUM, BIGNUM )
00226       SMLNUM = SQRT( SMLNUM ) / EPS
00227       BIGNUM = ONE / SMLNUM
00228 *
00229 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00230 *
00231       ANRM = CLANGE( 'M', N, N, A, LDA, DUM )
00232       SCALEA = .FALSE.
00233       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00234          SCALEA = .TRUE.
00235          CSCALE = SMLNUM
00236       ELSE IF( ANRM.GT.BIGNUM ) THEN
00237          SCALEA = .TRUE.
00238          CSCALE = BIGNUM
00239       END IF
00240       IF( SCALEA )
00241      $   CALL CLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
00242 *
00243 *     Permute the matrix to make it more nearly triangular
00244 *     (CWorkspace: none)
00245 *     (RWorkspace: need N)
00246 *
00247       IBAL = 1
00248       CALL CGEBAL( 'P', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
00249 *
00250 *     Reduce to upper Hessenberg form
00251 *     (CWorkspace: need 2*N, prefer N+N*NB)
00252 *     (RWorkspace: none)
00253 *
00254       ITAU = 1
00255       IWRK = N + ITAU
00256       CALL CGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
00257      $             LWORK-IWRK+1, IERR )
00258 *
00259       IF( WANTVS ) THEN
00260 *
00261 *        Copy Householder vectors to VS
00262 *
00263          CALL CLACPY( 'L', N, N, A, LDA, VS, LDVS )
00264 *
00265 *        Generate unitary matrix in VS
00266 *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
00267 *        (RWorkspace: none)
00268 *
00269          CALL CUNGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
00270      $                LWORK-IWRK+1, IERR )
00271       END IF
00272 *
00273       SDIM = 0
00274 *
00275 *     Perform QR iteration, accumulating Schur vectors in VS if desired
00276 *     (CWorkspace: need 1, prefer HSWORK (see comments) )
00277 *     (RWorkspace: none)
00278 *
00279       IWRK = ITAU
00280       CALL CHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, W, VS, LDVS,
00281      $             WORK( IWRK ), LWORK-IWRK+1, IEVAL )
00282       IF( IEVAL.GT.0 )
00283      $   INFO = IEVAL
00284 *
00285 *     Sort eigenvalues if desired
00286 *
00287       IF( WANTST .AND. INFO.EQ.0 ) THEN
00288          IF( SCALEA )
00289      $      CALL CLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, W, N, IERR )
00290          DO 10 I = 1, N
00291             BWORK( I ) = SELECT( W( I ) )
00292    10    CONTINUE
00293 *
00294 *        Reorder eigenvalues and transform Schur vectors
00295 *        (CWorkspace: none)
00296 *        (RWorkspace: none)
00297 *
00298          CALL CTRSEN( 'N', JOBVS, BWORK, N, A, LDA, VS, LDVS, W, SDIM,
00299      $                S, SEP, WORK( IWRK ), LWORK-IWRK+1, ICOND )
00300       END IF
00301 *
00302       IF( WANTVS ) THEN
00303 *
00304 *        Undo balancing
00305 *        (CWorkspace: none)
00306 *        (RWorkspace: need N)
00307 *
00308          CALL CGEBAK( 'P', 'R', N, ILO, IHI, RWORK( IBAL ), N, VS, LDVS,
00309      $                IERR )
00310       END IF
00311 *
00312       IF( SCALEA ) THEN
00313 *
00314 *        Undo scaling for the Schur form of A
00315 *
00316          CALL CLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
00317          CALL CCOPY( N, A, LDA+1, W, 1 )
00318       END IF
00319 *
00320       WORK( 1 ) = MAXWRK
00321       RETURN
00322 *
00323 *     End of CGEES
00324 *
00325       END
 All Files Functions