LAPACK: SRC/zgges.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
00002      $                  SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
00003      $                  LWORK, RWORK, BWORK, 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, SORT
00012       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
00013 *     ..
00014 *     .. Array Arguments ..
00015       LOGICAL            BWORK( * )
00016       DOUBLE PRECISION   RWORK( * )
00017       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
00018      $                   BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
00019      $                   WORK( * )
00020 *     ..
00021 *     .. Function Arguments ..
00022       LOGICAL            SELCTG
00023       EXTERNAL           SELCTG
00024 *     ..
00025 *
00026 *  Purpose
00027 *  =======
00028 *
00029 *  ZGGES computes for a pair of N-by-N complex nonsymmetric matrices
00030 *  (A,B), the generalized eigenvalues, the generalized complex Schur
00031 *  form (S, T), and optionally left and/or right Schur vectors (VSL
00032 *  and VSR). This gives the generalized Schur factorization
00033 *
00034 *          (A,B) = ( (VSL)*S*(VSR)**H, (VSL)*T*(VSR)**H )
00035 *
00036 *  where (VSR)**H is the conjugate-transpose of VSR.
00037 *
00038 *  Optionally, it also orders the eigenvalues so that a selected cluster
00039 *  of eigenvalues appears in the leading diagonal blocks of the upper
00040 *  triangular matrix S and the upper triangular matrix T. The leading
00041 *  columns of VSL and VSR then form an unitary basis for the
00042 *  corresponding left and right eigenspaces (deflating subspaces).
00043 *
00044 *  (If only the generalized eigenvalues are needed, use the driver
00045 *  ZGGEV instead, which is faster.)
00046 *
00047 *  A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
00048 *  or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
00049 *  usually represented as the pair (alpha,beta), as there is a
00050 *  reasonable interpretation for beta=0, and even for both being zero.
00051 *
00052 *  A pair of matrices (S,T) is in generalized complex Schur form if S
00053 *  and T are upper triangular and, in addition, the diagonal elements
00054 *  of T are non-negative real numbers.
00055 *
00056 *  Arguments
00057 *  =========
00058 *
00059 *  JOBVSL  (input) CHARACTER*1
00060 *          = 'N':  do not compute the left Schur vectors;
00061 *          = 'V':  compute the left Schur vectors.
00062 *
00063 *  JOBVSR  (input) CHARACTER*1
00064 *          = 'N':  do not compute the right Schur vectors;
00065 *          = 'V':  compute the right Schur vectors.
00066 *
00067 *  SORT    (input) CHARACTER*1
00068 *          Specifies whether or not to order the eigenvalues on the
00069 *          diagonal of the generalized Schur form.
00070 *          = 'N':  Eigenvalues are not ordered;
00071 *          = 'S':  Eigenvalues are ordered (see SELCTG).
00072 *
00073 *  SELCTG  (external procedure) LOGICAL FUNCTION of two COMPLEX*16 arguments
00074 *          SELCTG must be declared EXTERNAL in the calling subroutine.
00075 *          If SORT = 'N', SELCTG is not referenced.
00076 *          If SORT = 'S', SELCTG is used to select eigenvalues to sort
00077 *          to the top left of the Schur form.
00078 *          An eigenvalue ALPHA(j)/BETA(j) is selected if
00079 *          SELCTG(ALPHA(j),BETA(j)) is true.
00080 *
00081 *          Note that a selected complex eigenvalue may no longer satisfy
00082 *          SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since
00083 *          ordering may change the value of complex eigenvalues
00084 *          (especially if the eigenvalue is ill-conditioned), in this
00085 *          case INFO is set to N+2 (See INFO below).
00086 *
00087 *  N       (input) INTEGER
00088 *          The order of the matrices A, B, VSL, and VSR.  N >= 0.
00089 *
00090 *  A       (input/output) COMPLEX*16 array, dimension (LDA, N)
00091 *          On entry, the first of the pair of matrices.
00092 *          On exit, A has been overwritten by its generalized Schur
00093 *          form S.
00094 *
00095 *  LDA     (input) INTEGER
00096 *          The leading dimension of A.  LDA >= max(1,N).
00097 *
00098 *  B       (input/output) COMPLEX*16 array, dimension (LDB, N)
00099 *          On entry, the second of the pair of matrices.
00100 *          On exit, B has been overwritten by its generalized Schur
00101 *          form T.
00102 *
00103 *  LDB     (input) INTEGER
00104 *          The leading dimension of B.  LDB >= max(1,N).
00105 *
00106 *  SDIM    (output) INTEGER
00107 *          If SORT = 'N', SDIM = 0.
00108 *          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
00109 *          for which SELCTG is true.
00110 *
00111 *  ALPHA   (output) COMPLEX*16 array, dimension (N)
00112 *  BETA    (output) COMPLEX*16 array, dimension (N)
00113 *          On exit,  ALPHA(j)/BETA(j), j=1,...,N, will be the
00114 *          generalized eigenvalues.  ALPHA(j), j=1,...,N  and  BETA(j),
00115 *          j=1,...,N  are the diagonals of the complex Schur form (A,B)
00116 *          output by ZGGES. The  BETA(j) will be non-negative real.
00117 *
00118 *          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
00119 *          underflow, and BETA(j) may even be zero.  Thus, the user
00120 *          should avoid naively computing the ratio alpha/beta.
00121 *          However, ALPHA will be always less than and usually
00122 *          comparable with norm(A) in magnitude, and BETA always less
00123 *          than and usually comparable with norm(B).
00124 *
00125 *  VSL     (output) COMPLEX*16 array, dimension (LDVSL,N)
00126 *          If JOBVSL = 'V', VSL will contain the left Schur vectors.
00127 *          Not referenced if JOBVSL = 'N'.
00128 *
00129 *  LDVSL   (input) INTEGER
00130 *          The leading dimension of the matrix VSL. LDVSL >= 1, and
00131 *          if JOBVSL = 'V', LDVSL >= N.
00132 *
00133 *  VSR     (output) COMPLEX*16 array, dimension (LDVSR,N)
00134 *          If JOBVSR = 'V', VSR will contain the right Schur vectors.
00135 *          Not referenced if JOBVSR = 'N'.
00136 *
00137 *  LDVSR   (input) INTEGER
00138 *          The leading dimension of the matrix VSR. LDVSR >= 1, and
00139 *          if JOBVSR = 'V', LDVSR >= N.
00140 *
00141 *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
00142 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00143 *
00144 *  LWORK   (input) INTEGER
00145 *          The dimension of the array WORK.  LWORK >= max(1,2*N).
00146 *          For good performance, LWORK must generally be larger.
00147 *
00148 *          If LWORK = -1, then a workspace query is assumed; the routine
00149 *          only calculates the optimal size of the WORK array, returns
00150 *          this value as the first entry of the WORK array, and no error
00151 *          message related to LWORK is issued by XERBLA.
00152 *
00153 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (8*N)
00154 *
00155 *  BWORK   (workspace) LOGICAL array, dimension (N)
00156 *          Not referenced if SORT = 'N'.
00157 *
00158 *  INFO    (output) INTEGER
00159 *          = 0:  successful exit
00160 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00161 *          =1,...,N:
00162 *                The QZ iteration failed.  (A,B) are not in Schur
00163 *                form, but ALPHA(j) and BETA(j) should be correct for
00164 *                j=INFO+1,...,N.
00165 *          > N:  =N+1: other than QZ iteration failed in ZHGEQZ
00166 *                =N+2: after reordering, roundoff changed values of
00167 *                      some complex eigenvalues so that leading
00168 *                      eigenvalues in the Generalized Schur form no
00169 *                      longer satisfy SELCTG=.TRUE.  This could also
00170 *                      be caused due to scaling.
00171 *                =N+3: reordering falied in ZTGSEN.
00172 *
00173 *  =====================================================================
00174 *
00175 *     .. Parameters ..
00176       DOUBLE PRECISION   ZERO, ONE
00177       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00178       COMPLEX*16         CZERO, CONE
00179       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
00180      $                   CONE = ( 1.0D0, 0.0D0 ) )
00181 *     ..
00182 *     .. Local Scalars ..
00183       LOGICAL            CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
00184      $                   LQUERY, WANTST
00185       INTEGER            I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
00186      $                   ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, LWKMIN,
00187      $                   LWKOPT
00188       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
00189      $                   PVSR, SMLNUM
00190 *     ..
00191 *     .. Local Arrays ..
00192       INTEGER            IDUM( 1 )
00193       DOUBLE PRECISION   DIF( 2 )
00194 *     ..
00195 *     .. External Subroutines ..
00196       EXTERNAL           DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD,
00197      $                   ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGSEN, ZUNGQR,
00198      $                   ZUNMQR
00199 *     ..
00200 *     .. External Functions ..
00201       LOGICAL            LSAME
00202       INTEGER            ILAENV
00203       DOUBLE PRECISION   DLAMCH, ZLANGE
00204       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
00205 *     ..
00206 *     .. Intrinsic Functions ..
00207       INTRINSIC          MAX, SQRT
00208 *     ..
00209 *     .. Executable Statements ..
00210 *
00211 *     Decode the input arguments
00212 *
00213       IF( LSAME( JOBVSL, 'N' ) ) THEN
00214          IJOBVL = 1
00215          ILVSL = .FALSE.
00216       ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
00217          IJOBVL = 2
00218          ILVSL = .TRUE.
00219       ELSE
00220          IJOBVL = -1
00221          ILVSL = .FALSE.
00222       END IF
00223 *
00224       IF( LSAME( JOBVSR, 'N' ) ) THEN
00225          IJOBVR = 1
00226          ILVSR = .FALSE.
00227       ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
00228          IJOBVR = 2
00229          ILVSR = .TRUE.
00230       ELSE
00231          IJOBVR = -1
00232          ILVSR = .FALSE.
00233       END IF
00234 *
00235       WANTST = LSAME( SORT, 'S' )
00236 *
00237 *     Test the input arguments
00238 *
00239       INFO = 0
00240       LQUERY = ( LWORK.EQ.-1 )
00241       IF( IJOBVL.LE.0 ) THEN
00242          INFO = -1
00243       ELSE IF( IJOBVR.LE.0 ) THEN
00244          INFO = -2
00245       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
00246          INFO = -3
00247       ELSE IF( N.LT.0 ) THEN
00248          INFO = -5
00249       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00250          INFO = -7
00251       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00252          INFO = -9
00253       ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
00254          INFO = -14
00255       ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
00256          INFO = -16
00257       END IF
00258 *
00259 *     Compute workspace
00260 *      (Note: Comments in the code beginning "Workspace:" describe the
00261 *       minimal amount of workspace needed at that point in the code,
00262 *       as well as the preferred amount for good performance.
00263 *       NB refers to the optimal block size for the immediately
00264 *       following subroutine, as returned by ILAENV.)
00265 *
00266       IF( INFO.EQ.0 ) THEN
00267          LWKMIN = MAX( 1, 2*N )
00268          LWKOPT = MAX( 1, N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
00269          LWKOPT = MAX( LWKOPT, N +
00270      $                 N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, -1 ) )
00271          IF( ILVSL ) THEN
00272             LWKOPT = MAX( LWKOPT, N +
00273      $                    N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) )
00274          END IF
00275          WORK( 1 ) = LWKOPT
00276 *
00277          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY )
00278      $      INFO = -18
00279       END IF
00280 *
00281       IF( INFO.NE.0 ) THEN
00282          CALL XERBLA( 'ZGGES ', -INFO )
00283          RETURN
00284       ELSE IF( LQUERY ) THEN
00285          RETURN
00286       END IF
00287 *
00288 *     Quick return if possible
00289 *
00290       IF( N.EQ.0 ) THEN
00291          SDIM = 0
00292          RETURN
00293       END IF
00294 *
00295 *     Get machine constants
00296 *
00297       EPS = DLAMCH( 'P' )
00298       SMLNUM = DLAMCH( 'S' )
00299       BIGNUM = ONE / SMLNUM
00300       CALL DLABAD( SMLNUM, BIGNUM )
00301       SMLNUM = SQRT( SMLNUM ) / EPS
00302       BIGNUM = ONE / SMLNUM
00303 *
00304 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00305 *
00306       ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
00307       ILASCL = .FALSE.
00308       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00309          ANRMTO = SMLNUM
00310          ILASCL = .TRUE.
00311       ELSE IF( ANRM.GT.BIGNUM ) THEN
00312          ANRMTO = BIGNUM
00313          ILASCL = .TRUE.
00314       END IF
00315 *
00316       IF( ILASCL )
00317      $   CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
00318 *
00319 *     Scale B if max element outside range [SMLNUM,BIGNUM]
00320 *
00321       BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
00322       ILBSCL = .FALSE.
00323       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00324          BNRMTO = SMLNUM
00325          ILBSCL = .TRUE.
00326       ELSE IF( BNRM.GT.BIGNUM ) THEN
00327          BNRMTO = BIGNUM
00328          ILBSCL = .TRUE.
00329       END IF
00330 *
00331       IF( ILBSCL )
00332      $   CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
00333 *
00334 *     Permute the matrix to make it more nearly triangular
00335 *     (Real Workspace: need 6*N)
00336 *
00337       ILEFT = 1
00338       IRIGHT = N + 1
00339       IRWRK = IRIGHT + N
00340       CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
00341      $             RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
00342 *
00343 *     Reduce B to triangular form (QR decomposition of B)
00344 *     (Complex Workspace: need N, prefer N*NB)
00345 *
00346       IROWS = IHI + 1 - ILO
00347       ICOLS = N + 1 - ILO
00348       ITAU = 1
00349       IWRK = ITAU + IROWS
00350       CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
00351      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
00352 *
00353 *     Apply the orthogonal transformation to matrix A
00354 *     (Complex Workspace: need N, prefer N*NB)
00355 *
00356       CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
00357      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
00358      $             LWORK+1-IWRK, IERR )
00359 *
00360 *     Initialize VSL
00361 *     (Complex Workspace: need N, prefer N*NB)
00362 *
00363       IF( ILVSL ) THEN
00364          CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
00365          IF( IROWS.GT.1 ) THEN
00366             CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
00367      $                   VSL( ILO+1, ILO ), LDVSL )
00368          END IF
00369          CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
00370      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
00371       END IF
00372 *
00373 *     Initialize VSR
00374 *
00375       IF( ILVSR )
00376      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
00377 *
00378 *     Reduce to generalized Hessenberg form
00379 *     (Workspace: none needed)
00380 *
00381       CALL ZGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
00382      $             LDVSL, VSR, LDVSR, IERR )
00383 *
00384       SDIM = 0
00385 *
00386 *     Perform QZ algorithm, computing Schur vectors if desired
00387 *     (Complex Workspace: need N)
00388 *     (Real Workspace: need N)
00389 *
00390       IWRK = ITAU
00391       CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
00392      $             ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ),
00393      $             LWORK+1-IWRK, RWORK( IRWRK ), IERR )
00394       IF( IERR.NE.0 ) THEN
00395          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
00396             INFO = IERR
00397          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
00398             INFO = IERR - N
00399          ELSE
00400             INFO = N + 1
00401          END IF
00402          GO TO 30
00403       END IF
00404 *
00405 *     Sort eigenvalues ALPHA/BETA if desired
00406 *     (Workspace: none needed)
00407 *
00408       IF( WANTST ) THEN
00409 *
00410 *        Undo scaling on eigenvalues before selecting
00411 *
00412          IF( ILASCL )
00413      $      CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, 1, ALPHA, N, IERR )
00414          IF( ILBSCL )
00415      $      CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, 1, BETA, N, IERR )
00416 *
00417 *        Select eigenvalues
00418 *
00419          DO 10 I = 1, N
00420             BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) )
00421    10    CONTINUE
00422 *
00423          CALL ZTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHA,
00424      $                BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL, PVSR,
00425      $                DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1, IERR )
00426          IF( IERR.EQ.1 )
00427      $      INFO = N + 3
00428 *
00429       END IF
00430 *
00431 *     Apply back-permutation to VSL and VSR
00432 *     (Workspace: none needed)
00433 *
00434       IF( ILVSL )
00435      $   CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
00436      $                RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
00437       IF( ILVSR )
00438      $   CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
00439      $                RWORK( IRIGHT ), N, VSR, LDVSR, IERR )
00440 *
00441 *     Undo scaling
00442 *
00443       IF( ILASCL ) THEN
00444          CALL ZLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
00445          CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
00446       END IF
00447 *
00448       IF( ILBSCL ) THEN
00449          CALL ZLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
00450          CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
00451       END IF
00452 *
00453       IF( WANTST ) THEN
00454 *
00455 *        Check if reordering is correct
00456 *
00457          LASTSL = .TRUE.
00458          SDIM = 0
00459          DO 20 I = 1, N
00460             CURSL = SELCTG( ALPHA( I ), BETA( I ) )
00461             IF( CURSL )
00462      $         SDIM = SDIM + 1
00463             IF( CURSL .AND. .NOT.LASTSL )
00464      $         INFO = N + 2
00465             LASTSL = CURSL
00466    20    CONTINUE
00467 *
00468       END IF
00469 *
00470    30 CONTINUE
00471 *
00472       WORK( 1 ) = LWKOPT
00473 *
00474       RETURN
00475 *
00476 *     End of ZGGES
00477 *
00478       END