LAPACK: SRC/zgegs.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZGEGS( 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       DOUBLE PRECISION   RWORK( * )
00016       COMPLEX*16         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 ZGGES.
00025 *
00026 *  ZGEGS 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 *  ZGEGV should be used instead.  See ZGEGV 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*16 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*16 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*16 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*16 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 *
00089 *  VSL     (output) COMPLEX*16 array, dimension (LDVSL,N)
00090 *          If JOBVSL = 'V', the matrix of left Schur vectors Q.
00091 *          Not referenced if JOBVSL = 'N'.
00092 *
00093 *  LDVSL   (input) INTEGER
00094 *          The leading dimension of the matrix VSL. LDVSL >= 1, and
00095 *          if JOBVSL = 'V', LDVSL >= N.
00096 *
00097 *  VSR     (output) COMPLEX*16 array, dimension (LDVSR,N)
00098 *          If JOBVSR = 'V', the matrix of right Schur vectors Z.
00099 *          Not referenced if JOBVSR = 'N'.
00100 *
00101 *  LDVSR   (input) INTEGER
00102 *          The leading dimension of the matrix VSR. LDVSR >= 1, and
00103 *          if JOBVSR = 'V', LDVSR >= N.
00104 *
00105 *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
00106 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00107 *
00108 *  LWORK   (input) INTEGER
00109 *          The dimension of the array WORK.  LWORK >= max(1,2*N).
00110 *          For good performance, LWORK must generally be larger.
00111 *          To compute the optimal value of LWORK, call ILAENV to get
00112 *          blocksizes (for ZGEQRF, ZUNMQR, and CUNGQR.)  Then compute:
00113 *          NB  -- MAX of the blocksizes for ZGEQRF, ZUNMQR, and CUNGQR;
00114 *          the optimal LWORK is N*(NB+1).
00115 *
00116 *          If LWORK = -1, then a workspace query is assumed; the routine
00117 *          only calculates the optimal size of the WORK array, returns
00118 *          this value as the first entry of the WORK array, and no error
00119 *          message related to LWORK is issued by XERBLA.
00120 *
00121 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (3*N)
00122 *
00123 *  INFO    (output) INTEGER
00124 *          = 0:  successful exit
00125 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00126 *          =1,...,N:
00127 *                The QZ iteration failed.  (A,B) are not in Schur
00128 *                form, but ALPHA(j) and BETA(j) should be correct for
00129 *                j=INFO+1,...,N.
00130 *          > N:  errors that usually indicate LAPACK problems:
00131 *                =N+1: error return from ZGGBAL
00132 *                =N+2: error return from ZGEQRF
00133 *                =N+3: error return from ZUNMQR
00134 *                =N+4: error return from ZUNGQR
00135 *                =N+5: error return from ZGGHRD
00136 *                =N+6: error return from ZHGEQZ (other than failed
00137 *                                               iteration)
00138 *                =N+7: error return from ZGGBAK (computing VSL)
00139 *                =N+8: error return from ZGGBAK (computing VSR)
00140 *                =N+9: error return from ZLASCL (various places)
00141 *
00142 *  =====================================================================
00143 *
00144 *     .. Parameters ..
00145       DOUBLE PRECISION   ZERO, ONE
00146       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00147       COMPLEX*16         CZERO, CONE
00148       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
00149      $                   CONE = ( 1.0D0, 0.0D0 ) )
00150 *     ..
00151 *     .. Local Scalars ..
00152       LOGICAL            ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
00153       INTEGER            ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
00154      $                   IRIGHT, IROWS, IRWORK, ITAU, IWORK, LOPT,
00155      $                   LWKMIN, LWKOPT, NB, NB1, NB2, NB3
00156       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
00157      $                   SAFMIN, SMLNUM
00158 *     ..
00159 *     .. External Subroutines ..
00160       EXTERNAL           XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD, ZHGEQZ,
00161      $                   ZLACPY, ZLASCL, ZLASET, ZUNGQR, ZUNMQR
00162 *     ..
00163 *     .. External Functions ..
00164       LOGICAL            LSAME
00165       INTEGER            ILAENV
00166       DOUBLE PRECISION   DLAMCH, ZLANGE
00167       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
00168 *     ..
00169 *     .. Intrinsic Functions ..
00170       INTRINSIC          INT, MAX
00171 *     ..
00172 *     .. Executable Statements ..
00173 *
00174 *     Decode the input arguments
00175 *
00176       IF( LSAME( JOBVSL, 'N' ) ) THEN
00177          IJOBVL = 1
00178          ILVSL = .FALSE.
00179       ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
00180          IJOBVL = 2
00181          ILVSL = .TRUE.
00182       ELSE
00183          IJOBVL = -1
00184          ILVSL = .FALSE.
00185       END IF
00186 *
00187       IF( LSAME( JOBVSR, 'N' ) ) THEN
00188          IJOBVR = 1
00189          ILVSR = .FALSE.
00190       ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
00191          IJOBVR = 2
00192          ILVSR = .TRUE.
00193       ELSE
00194          IJOBVR = -1
00195          ILVSR = .FALSE.
00196       END IF
00197 *
00198 *     Test the input arguments
00199 *
00200       LWKMIN = MAX( 2*N, 1 )
00201       LWKOPT = LWKMIN
00202       WORK( 1 ) = LWKOPT
00203       LQUERY = ( LWORK.EQ.-1 )
00204       INFO = 0
00205       IF( IJOBVL.LE.0 ) THEN
00206          INFO = -1
00207       ELSE IF( IJOBVR.LE.0 ) THEN
00208          INFO = -2
00209       ELSE IF( N.LT.0 ) THEN
00210          INFO = -3
00211       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00212          INFO = -5
00213       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00214          INFO = -7
00215       ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
00216          INFO = -11
00217       ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
00218          INFO = -13
00219       ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
00220          INFO = -15
00221       END IF
00222 *
00223       IF( INFO.EQ.0 ) THEN
00224          NB1 = ILAENV( 1, 'ZGEQRF', ' ', N, N, -1, -1 )
00225          NB2 = ILAENV( 1, 'ZUNMQR', ' ', N, N, N, -1 )
00226          NB3 = ILAENV( 1, 'ZUNGQR', ' ', N, N, N, -1 )
00227          NB = MAX( NB1, NB2, NB3 )
00228          LOPT = N*( NB+1 )
00229          WORK( 1 ) = LOPT
00230       END IF
00231 *
00232       IF( INFO.NE.0 ) THEN
00233          CALL XERBLA( 'ZGEGS ', -INFO )
00234          RETURN
00235       ELSE IF( LQUERY ) THEN
00236          RETURN
00237       END IF
00238 *
00239 *     Quick return if possible
00240 *
00241       IF( N.EQ.0 )
00242      $   RETURN
00243 *
00244 *     Get machine constants
00245 *
00246       EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
00247       SAFMIN = DLAMCH( 'S' )
00248       SMLNUM = N*SAFMIN / EPS
00249       BIGNUM = ONE / SMLNUM
00250 *
00251 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00252 *
00253       ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
00254       ILASCL = .FALSE.
00255       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00256          ANRMTO = SMLNUM
00257          ILASCL = .TRUE.
00258       ELSE IF( ANRM.GT.BIGNUM ) THEN
00259          ANRMTO = BIGNUM
00260          ILASCL = .TRUE.
00261       END IF
00262 *
00263       IF( ILASCL ) THEN
00264          CALL ZLASCL( 'G', -1, -1, ANRM, ANRMTO, N, N, A, LDA, IINFO )
00265          IF( IINFO.NE.0 ) THEN
00266             INFO = N + 9
00267             RETURN
00268          END IF
00269       END IF
00270 *
00271 *     Scale B if max element outside range [SMLNUM,BIGNUM]
00272 *
00273       BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
00274       ILBSCL = .FALSE.
00275       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00276          BNRMTO = SMLNUM
00277          ILBSCL = .TRUE.
00278       ELSE IF( BNRM.GT.BIGNUM ) THEN
00279          BNRMTO = BIGNUM
00280          ILBSCL = .TRUE.
00281       END IF
00282 *
00283       IF( ILBSCL ) THEN
00284          CALL ZLASCL( 'G', -1, -1, BNRM, BNRMTO, N, N, B, LDB, IINFO )
00285          IF( IINFO.NE.0 ) THEN
00286             INFO = N + 9
00287             RETURN
00288          END IF
00289       END IF
00290 *
00291 *     Permute the matrix to make it more nearly triangular
00292 *
00293       ILEFT = 1
00294       IRIGHT = N + 1
00295       IRWORK = IRIGHT + N
00296       IWORK = 1
00297       CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
00298      $             RWORK( IRIGHT ), RWORK( IRWORK ), IINFO )
00299       IF( IINFO.NE.0 ) THEN
00300          INFO = N + 1
00301          GO TO 10
00302       END IF
00303 *
00304 *     Reduce B to triangular form, and initialize VSL and/or VSR
00305 *
00306       IROWS = IHI + 1 - ILO
00307       ICOLS = N + 1 - ILO
00308       ITAU = IWORK
00309       IWORK = ITAU + IROWS
00310       CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
00311      $             WORK( IWORK ), LWORK+1-IWORK, IINFO )
00312       IF( IINFO.GE.0 )
00313      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00314       IF( IINFO.NE.0 ) THEN
00315          INFO = N + 2
00316          GO TO 10
00317       END IF
00318 *
00319       CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
00320      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
00321      $             LWORK+1-IWORK, IINFO )
00322       IF( IINFO.GE.0 )
00323      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00324       IF( IINFO.NE.0 ) THEN
00325          INFO = N + 3
00326          GO TO 10
00327       END IF
00328 *
00329       IF( ILVSL ) THEN
00330          CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
00331          CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
00332      $                VSL( ILO+1, ILO ), LDVSL )
00333          CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
00334      $                WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
00335      $                IINFO )
00336          IF( IINFO.GE.0 )
00337      $      LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00338          IF( IINFO.NE.0 ) THEN
00339             INFO = N + 4
00340             GO TO 10
00341          END IF
00342       END IF
00343 *
00344       IF( ILVSR )
00345      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
00346 *
00347 *     Reduce to generalized Hessenberg form
00348 *
00349       CALL ZGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
00350      $             LDVSL, VSR, LDVSR, IINFO )
00351       IF( IINFO.NE.0 ) THEN
00352          INFO = N + 5
00353          GO TO 10
00354       END IF
00355 *
00356 *     Perform QZ algorithm, computing Schur vectors if desired
00357 *
00358       IWORK = ITAU
00359       CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
00360      $             ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWORK ),
00361      $             LWORK+1-IWORK, RWORK( IRWORK ), IINFO )
00362       IF( IINFO.GE.0 )
00363      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00364       IF( IINFO.NE.0 ) THEN
00365          IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
00366             INFO = IINFO
00367          ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
00368             INFO = IINFO - N
00369          ELSE
00370             INFO = N + 6
00371          END IF
00372          GO TO 10
00373       END IF
00374 *
00375 *     Apply permutation to VSL and VSR
00376 *
00377       IF( ILVSL ) THEN
00378          CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
00379      $                RWORK( IRIGHT ), N, VSL, LDVSL, IINFO )
00380          IF( IINFO.NE.0 ) THEN
00381             INFO = N + 7
00382             GO TO 10
00383          END IF
00384       END IF
00385       IF( ILVSR ) THEN
00386          CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
00387      $                RWORK( IRIGHT ), N, VSR, LDVSR, IINFO )
00388          IF( IINFO.NE.0 ) THEN
00389             INFO = N + 8
00390             GO TO 10
00391          END IF
00392       END IF
00393 *
00394 *     Undo scaling
00395 *
00396       IF( ILASCL ) THEN
00397          CALL ZLASCL( 'U', -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
00398          IF( IINFO.NE.0 ) THEN
00399             INFO = N + 9
00400             RETURN
00401          END IF
00402          CALL ZLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHA, N, IINFO )
00403          IF( IINFO.NE.0 ) THEN
00404             INFO = N + 9
00405             RETURN
00406          END IF
00407       END IF
00408 *
00409       IF( ILBSCL ) THEN
00410          CALL ZLASCL( 'U', -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
00411          IF( IINFO.NE.0 ) THEN
00412             INFO = N + 9
00413             RETURN
00414          END IF
00415          CALL ZLASCL( 'G', -1, -1, BNRMTO, BNRM, N, 1, BETA, N, IINFO )
00416          IF( IINFO.NE.0 ) THEN
00417             INFO = N + 9
00418             RETURN
00419          END IF
00420       END IF
00421 *
00422    10 CONTINUE
00423       WORK( 1 ) = LWKOPT
00424 *
00425       RETURN
00426 *
00427 *     End of ZGEGS
00428 *
00429       END