LAPACK: SRC/cgegs.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHA, BETA,
00002      $                  VSL, LDVSL, VSR, LDVSR, WORK, LWORK, RWORK,
00003      $                  INFO )
00004 *
00005 *  -- LAPACK driver routine (version 3.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     November 2006
00009 *
00010 *     .. Scalar Arguments ..
00011       CHARACTER          JOBVSL, JOBVSR
00012       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
00013 *     ..
00014 *     .. Array Arguments ..
00015       REAL               RWORK( * )
00016       COMPLEX            A( LDA, * ), ALPHA( * ), B( LDB, * ),
00017      $                   BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
00018      $                   WORK( * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  This routine is deprecated and has been replaced by routine CGGES.
00025 *
00026 *  CGEGS computes the eigenvalues, Schur form, and, optionally, the
00027 *  left and or/right Schur vectors of a complex matrix pair (A,B).
00028 *  Given two square matrices A and B, the generalized Schur
00029 *  factorization has the form
00030 *  
00031 *     A = Q*S*Z**H,  B = Q*T*Z**H
00032 *  
00033 *  where Q and Z are unitary matrices and S and T are upper triangular.
00034 *  The columns of Q are the left Schur vectors
00035 *  and the columns of Z are the right Schur vectors.
00036 *  
00037 *  If only the eigenvalues of (A,B) are needed, the driver routine
00038 *  CGEGV should be used instead.  See CGEGV for a description of the
00039 *  eigenvalues of the generalized nonsymmetric eigenvalue problem
00040 *  (GNEP).
00041 *
00042 *  Arguments
00043 *  =========
00044 *
00045 *  JOBVSL   (input) CHARACTER*1
00046 *          = 'N':  do not compute the left Schur vectors;
00047 *          = 'V':  compute the left Schur vectors (returned in VSL).
00048 *
00049 *  JOBVSR   (input) CHARACTER*1
00050 *          = 'N':  do not compute the right Schur vectors;
00051 *          = 'V':  compute the right Schur vectors (returned in VSR).
00052 *
00053 *  N       (input) INTEGER
00054 *          The order of the matrices A, B, VSL, and VSR.  N >= 0.
00055 *
00056 *  A       (input/output) COMPLEX array, dimension (LDA, N)
00057 *          On entry, the matrix A.
00058 *          On exit, the upper triangular matrix S from the generalized
00059 *          Schur factorization.
00060 *
00061 *  LDA     (input) INTEGER
00062 *          The leading dimension of A.  LDA >= max(1,N).
00063 *
00064 *  B       (input/output) COMPLEX array, dimension (LDB, N)
00065 *          On entry, the matrix B.
00066 *          On exit, the upper triangular matrix T from the generalized
00067 *          Schur factorization.
00068 *
00069 *  LDB     (input) INTEGER
00070 *          The leading dimension of B.  LDB >= max(1,N).
00071 *
00072 *  ALPHA   (output) COMPLEX array, dimension (N)
00073 *          The complex scalars alpha that define the eigenvalues of
00074 *          GNEP.  ALPHA(j) = S(j,j), the diagonal element of the Schur
00075 *          form of A.
00076 *
00077 *  BETA    (output) COMPLEX array, dimension (N)
00078 *          The non-negative real scalars beta that define the
00079 *          eigenvalues of GNEP.  BETA(j) = T(j,j), the diagonal element
00080 *          of the triangular factor T.
00081 *
00082 *          Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
00083 *          represent the j-th eigenvalue of the matrix pair (A,B), in
00084 *          one of the forms lambda = alpha/beta or mu = beta/alpha.
00085 *          Since either lambda or mu may overflow, they should not,
00086 *          in general, be computed.
00087 *
00088 *  VSL     (output) COMPLEX array, dimension (LDVSL,N)
00089 *          If JOBVSL = 'V', the matrix of left Schur vectors Q.
00090 *          Not referenced if JOBVSL = 'N'.
00091 *
00092 *  LDVSL   (input) INTEGER
00093 *          The leading dimension of the matrix VSL. LDVSL >= 1, and
00094 *          if JOBVSL = 'V', LDVSL >= N.
00095 *
00096 *  VSR     (output) COMPLEX array, dimension (LDVSR,N)
00097 *          If JOBVSR = 'V', the matrix of right Schur vectors Z.
00098 *          Not referenced if JOBVSR = 'N'.
00099 *
00100 *  LDVSR   (input) INTEGER
00101 *          The leading dimension of the matrix VSR. LDVSR >= 1, and
00102 *          if JOBVSR = 'V', LDVSR >= N.
00103 *
00104 *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
00105 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00106 *
00107 *  LWORK   (input) INTEGER
00108 *          The dimension of the array WORK.  LWORK >= max(1,2*N).
00109 *          For good performance, LWORK must generally be larger.
00110 *          To compute the optimal value of LWORK, call ILAENV to get
00111 *          blocksizes (for CGEQRF, CUNMQR, and CUNGQR.)  Then compute:
00112 *          NB  -- MAX of the blocksizes for CGEQRF, CUNMQR, and CUNGQR;
00113 *          the optimal LWORK is N*(NB+1).
00114 *
00115 *          If LWORK = -1, then a workspace query is assumed; the routine
00116 *          only calculates the optimal size of the WORK array, returns
00117 *          this value as the first entry of the WORK array, and no error
00118 *          message related to LWORK is issued by XERBLA.
00119 *
00120 *  RWORK   (workspace) REAL array, dimension (3*N)
00121 *
00122 *  INFO    (output) INTEGER
00123 *          = 0:  successful exit
00124 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00125 *          =1,...,N:
00126 *                The QZ iteration failed.  (A,B) are not in Schur
00127 *                form, but ALPHA(j) and BETA(j) should be correct for
00128 *                j=INFO+1,...,N.
00129 *          > N:  errors that usually indicate LAPACK problems:
00130 *                =N+1: error return from CGGBAL
00131 *                =N+2: error return from CGEQRF
00132 *                =N+3: error return from CUNMQR
00133 *                =N+4: error return from CUNGQR
00134 *                =N+5: error return from CGGHRD
00135 *                =N+6: error return from CHGEQZ (other than failed
00136 *                                               iteration)
00137 *                =N+7: error return from CGGBAK (computing VSL)
00138 *                =N+8: error return from CGGBAK (computing VSR)
00139 *                =N+9: error return from CLASCL (various places)
00140 *
00141 *  =====================================================================
00142 *
00143 *     .. Parameters ..
00144       REAL               ZERO, ONE
00145       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00146       COMPLEX            CZERO, CONE
00147       PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ),
00148      $                   CONE = ( 1.0E0, 0.0E0 ) )
00149 *     ..
00150 *     .. Local Scalars ..
00151       LOGICAL            ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
00152       INTEGER            ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT,
00153      $                   ILO, IRIGHT, IROWS, IRWORK, ITAU, IWORK,
00154      $                   LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3
00155       REAL               ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
00156      $                   SAFMIN, SMLNUM
00157 *     ..
00158 *     .. External Subroutines ..
00159       EXTERNAL           CGEQRF, CGGBAK, CGGBAL, CGGHRD, CHGEQZ, CLACPY,
00160      $                   CLASCL, CLASET, CUNGQR, CUNMQR, XERBLA
00161 *     ..
00162 *     .. External Functions ..
00163       LOGICAL            LSAME
00164       INTEGER            ILAENV
00165       REAL               CLANGE, SLAMCH
00166       EXTERNAL           ILAENV, LSAME, CLANGE, SLAMCH
00167 *     ..
00168 *     .. Intrinsic Functions ..
00169       INTRINSIC          INT, MAX
00170 *     ..
00171 *     .. Executable Statements ..
00172 *
00173 *     Decode the input arguments
00174 *
00175       IF( LSAME( JOBVSL, 'N' ) ) THEN
00176          IJOBVL = 1
00177          ILVSL = .FALSE.
00178       ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
00179          IJOBVL = 2
00180          ILVSL = .TRUE.
00181       ELSE
00182          IJOBVL = -1
00183          ILVSL = .FALSE.
00184       END IF
00185 *
00186       IF( LSAME( JOBVSR, 'N' ) ) THEN
00187          IJOBVR = 1
00188          ILVSR = .FALSE.
00189       ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
00190          IJOBVR = 2
00191          ILVSR = .TRUE.
00192       ELSE
00193          IJOBVR = -1
00194          ILVSR = .FALSE.
00195       END IF
00196 *
00197 *     Test the input arguments
00198 *
00199       LWKMIN = MAX( 2*N, 1 )
00200       LWKOPT = LWKMIN
00201       WORK( 1 ) = LWKOPT
00202       LQUERY = ( LWORK.EQ.-1 )
00203       INFO = 0
00204       IF( IJOBVL.LE.0 ) THEN
00205          INFO = -1
00206       ELSE IF( IJOBVR.LE.0 ) THEN
00207          INFO = -2
00208       ELSE IF( N.LT.0 ) THEN
00209          INFO = -3
00210       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00211          INFO = -5
00212       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00213          INFO = -7
00214       ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
00215          INFO = -11
00216       ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
00217          INFO = -13
00218       ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
00219          INFO = -15
00220       END IF
00221 *
00222       IF( INFO.EQ.0 ) THEN
00223          NB1 = ILAENV( 1, 'CGEQRF', ' ', N, N, -1, -1 )
00224          NB2 = ILAENV( 1, 'CUNMQR', ' ', N, N, N, -1 )
00225          NB3 = ILAENV( 1, 'CUNGQR', ' ', N, N, N, -1 )
00226          NB = MAX( NB1, NB2, NB3 )
00227          LOPT = N*(NB+1)
00228          WORK( 1 ) = LOPT
00229       END IF
00230 *
00231       IF( INFO.NE.0 ) THEN
00232          CALL XERBLA( 'CGEGS ', -INFO )
00233          RETURN
00234       ELSE IF( LQUERY ) THEN
00235          RETURN
00236       END IF
00237 *
00238 *     Quick return if possible
00239 *
00240       IF( N.EQ.0 )
00241      $   RETURN
00242 *
00243 *     Get machine constants
00244 *
00245       EPS = SLAMCH( 'E' )*SLAMCH( 'B' )
00246       SAFMIN = SLAMCH( 'S' )
00247       SMLNUM = N*SAFMIN / EPS
00248       BIGNUM = ONE / SMLNUM
00249 *
00250 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00251 *
00252       ANRM = CLANGE( 'M', N, N, A, LDA, RWORK )
00253       ILASCL = .FALSE.
00254       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00255          ANRMTO = SMLNUM
00256          ILASCL = .TRUE.
00257       ELSE IF( ANRM.GT.BIGNUM ) THEN
00258          ANRMTO = BIGNUM
00259          ILASCL = .TRUE.
00260       END IF
00261 *
00262       IF( ILASCL ) THEN
00263          CALL CLASCL( 'G', -1, -1, ANRM, ANRMTO, N, N, A, LDA, IINFO )
00264          IF( IINFO.NE.0 ) THEN
00265             INFO = N + 9
00266             RETURN
00267          END IF
00268       END IF
00269 *
00270 *     Scale B if max element outside range [SMLNUM,BIGNUM]
00271 *
00272       BNRM = CLANGE( 'M', N, N, B, LDB, RWORK )
00273       ILBSCL = .FALSE.
00274       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00275          BNRMTO = SMLNUM
00276          ILBSCL = .TRUE.
00277       ELSE IF( BNRM.GT.BIGNUM ) THEN
00278          BNRMTO = BIGNUM
00279          ILBSCL = .TRUE.
00280       END IF
00281 *
00282       IF( ILBSCL ) THEN
00283          CALL CLASCL( 'G', -1, -1, BNRM, BNRMTO, N, N, B, LDB, IINFO )
00284          IF( IINFO.NE.0 ) THEN
00285             INFO = N + 9
00286             RETURN
00287          END IF
00288       END IF
00289 *
00290 *     Permute the matrix to make it more nearly triangular
00291 *
00292       ILEFT = 1
00293       IRIGHT = N + 1
00294       IRWORK = IRIGHT + N
00295       IWORK = 1
00296       CALL CGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
00297      $             RWORK( IRIGHT ), RWORK( IRWORK ), IINFO )
00298       IF( IINFO.NE.0 ) THEN
00299          INFO = N + 1
00300          GO TO 10
00301       END IF
00302 *
00303 *     Reduce B to triangular form, and initialize VSL and/or VSR
00304 *
00305       IROWS = IHI + 1 - ILO
00306       ICOLS = N + 1 - ILO
00307       ITAU = IWORK
00308       IWORK = ITAU + IROWS
00309       CALL CGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
00310      $             WORK( IWORK ), LWORK+1-IWORK, IINFO )
00311       IF( IINFO.GE.0 )
00312      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00313       IF( IINFO.NE.0 ) THEN
00314          INFO = N + 2
00315          GO TO 10
00316       END IF
00317 *
00318       CALL CUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
00319      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
00320      $             LWORK+1-IWORK, IINFO )
00321       IF( IINFO.GE.0 )
00322      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00323       IF( IINFO.NE.0 ) THEN
00324          INFO = N + 3
00325          GO TO 10
00326       END IF
00327 *
00328       IF( ILVSL ) THEN
00329          CALL CLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
00330          CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
00331      $                VSL( ILO+1, ILO ), LDVSL )
00332          CALL CUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
00333      $                WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
00334      $                IINFO )
00335          IF( IINFO.GE.0 )
00336      $      LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00337          IF( IINFO.NE.0 ) THEN
00338             INFO = N + 4
00339             GO TO 10
00340          END IF
00341       END IF
00342 *
00343       IF( ILVSR )
00344      $   CALL CLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
00345 *
00346 *     Reduce to generalized Hessenberg form
00347 *
00348       CALL CGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
00349      $             LDVSL, VSR, LDVSR, IINFO )
00350       IF( IINFO.NE.0 ) THEN
00351          INFO = N + 5
00352          GO TO 10
00353       END IF
00354 *
00355 *     Perform QZ algorithm, computing Schur vectors if desired
00356 *
00357       IWORK = ITAU
00358       CALL CHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
00359      $             ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWORK ),
00360      $             LWORK+1-IWORK, RWORK( IRWORK ), IINFO )
00361       IF( IINFO.GE.0 )
00362      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00363       IF( IINFO.NE.0 ) THEN
00364          IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
00365             INFO = IINFO
00366          ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
00367             INFO = IINFO - N
00368          ELSE
00369             INFO = N + 6
00370          END IF
00371          GO TO 10
00372       END IF
00373 *
00374 *     Apply permutation to VSL and VSR
00375 *
00376       IF( ILVSL ) THEN
00377          CALL CGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
00378      $                RWORK( IRIGHT ), N, VSL, LDVSL, IINFO )
00379          IF( IINFO.NE.0 ) THEN
00380             INFO = N + 7
00381             GO TO 10
00382          END IF
00383       END IF
00384       IF( ILVSR ) THEN
00385          CALL CGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
00386      $                RWORK( IRIGHT ), N, VSR, LDVSR, IINFO )
00387          IF( IINFO.NE.0 ) THEN
00388             INFO = N + 8
00389             GO TO 10
00390          END IF
00391       END IF
00392 *
00393 *     Undo scaling
00394 *
00395       IF( ILASCL ) THEN
00396          CALL CLASCL( 'U', -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
00397          IF( IINFO.NE.0 ) THEN
00398             INFO = N + 9
00399             RETURN
00400          END IF
00401          CALL CLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHA, N, IINFO )
00402          IF( IINFO.NE.0 ) THEN
00403             INFO = N + 9
00404             RETURN
00405          END IF
00406       END IF
00407 *
00408       IF( ILBSCL ) THEN
00409          CALL CLASCL( 'U', -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
00410          IF( IINFO.NE.0 ) THEN
00411             INFO = N + 9
00412             RETURN
00413          END IF
00414          CALL CLASCL( 'G', -1, -1, BNRMTO, BNRM, N, 1, BETA, N, IINFO )
00415          IF( IINFO.NE.0 ) THEN
00416             INFO = N + 9
00417             RETURN
00418          END IF
00419       END IF
00420 *
00421    10 CONTINUE
00422       WORK( 1 ) = LWKOPT
00423 *
00424       RETURN
00425 *
00426 *     End of CGEGS
00427 *
00428       END