283 SUBROUTINE dgges( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
284 $ sdim, alphar, alphai, beta, vsl, ldvsl, vsr,
285 $ ldvsr, work, lwork, bwork, info )
293 CHARACTER JOBVSL, JOBVSR, SORT
294 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
298 DOUBLE PRECISION A( lda, * ), ALPHAI( * ), ALPHAR( * ),
299 $ b( ldb, * ), beta( * ), vsl( ldvsl, * ),
300 $ vsr( ldvsr, * ), work( * )
310 DOUBLE PRECISION ZERO, ONE
311 parameter( zero = 0.0d+0, one = 1.0d+0 )
314 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
315 $ lquery, lst2sl, wantst
316 INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
317 $ ilo, ip, iright, irows, itau, iwrk, maxwrk,
319 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
320 $ pvsr, safmax, safmin, smlnum
324 DOUBLE PRECISION DIF( 2 )
334 DOUBLE PRECISION DLAMCH, DLANGE
335 EXTERNAL lsame, ilaenv, dlamch, dlange
338 INTRINSIC abs, max, sqrt
344 IF( lsame( jobvsl,
'N' ) )
THEN
347 ELSE IF( lsame( jobvsl,
'V' ) )
THEN
355 IF( lsame( jobvsr,
'N' ) )
THEN
358 ELSE IF( lsame( jobvsr,
'V' ) )
THEN
366 wantst = lsame( sort,
'S' )
371 lquery = ( lwork.EQ.-1 )
372 IF( ijobvl.LE.0 )
THEN
374 ELSE IF( ijobvr.LE.0 )
THEN
376 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort,
'N' ) ) )
THEN
378 ELSE IF( n.LT.0 )
THEN
380 ELSE IF( lda.LT.max( 1, n ) )
THEN
382 ELSE IF( ldb.LT.max( 1, n ) )
THEN
384 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) )
THEN
386 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) )
THEN
399 minwrk = max( 8*n, 6*n + 16 )
400 maxwrk = minwrk - n +
401 $ n*ilaenv( 1,
'DGEQRF',
' ', n, 1, n, 0 )
402 maxwrk = max( maxwrk, minwrk - n +
403 $ n*ilaenv( 1,
'DORMQR',
' ', n, 1, n, -1 ) )
405 maxwrk = max( maxwrk, minwrk - n +
406 $ n*ilaenv( 1,
'DORGQR',
' ', n, 1, n, -1 ) )
414 IF( lwork.LT.minwrk .AND. .NOT.lquery )
419 CALL
xerbla(
'DGGES ', -info )
421 ELSE IF( lquery )
THEN
435 safmin = dlamch(
'S' )
436 safmax = one / safmin
437 CALL
dlabad( safmin, safmax )
438 smlnum = sqrt( safmin ) / eps
439 bignum = one / smlnum
443 anrm = dlange(
'M', n, n, a, lda, work )
445 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
448 ELSE IF( anrm.GT.bignum )
THEN
453 $ CALL
dlascl(
'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
457 bnrm = dlange(
'M', n, n, b, ldb, work )
459 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
462 ELSE IF( bnrm.GT.bignum )
THEN
467 $ CALL
dlascl(
'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
475 CALL
dggbal(
'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
476 $ work( iright ), work( iwrk ), ierr )
481 irows = ihi + 1 - ilo
485 CALL
dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
486 $ work( iwrk ), lwork+1-iwrk, ierr )
491 CALL
dormqr(
'L',
'T', irows, icols, irows, b( ilo, ilo ), ldb,
492 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
493 $ lwork+1-iwrk, ierr )
499 CALL
dlaset(
'Full', n, n, zero, one, vsl, ldvsl )
500 IF( irows.GT.1 )
THEN
501 CALL
dlacpy(
'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
502 $ vsl( ilo+1, ilo ), ldvsl )
504 CALL
dorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
505 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
511 $ CALL
dlaset(
'Full', n, n, zero, one, vsr, ldvsr )
516 CALL
dgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
517 $ ldvsl, vsr, ldvsr, ierr )
523 CALL
dhgeqz(
'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
524 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
525 $ work( iwrk ), lwork+1-iwrk, ierr )
527 IF( ierr.GT.0 .AND. ierr.LE.n )
THEN
529 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n )
THEN
546 CALL
dlascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
548 CALL
dlascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
552 $ CALL
dlascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
557 bwork( i ) = selctg( alphar( i ), alphai( i ), beta( i ) )
560 CALL
dtgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb, alphar,
561 $ alphai, beta, vsl, ldvsl, vsr, ldvsr, sdim, pvsl,
562 $ pvsr, dif, work( iwrk ), lwork-iwrk+1, idum, 1,
573 $ CALL
dggbak(
'P',
'L', n, ilo, ihi, work( ileft ),
574 $ work( iright ), n, vsl, ldvsl, ierr )
577 $ CALL
dggbak(
'P',
'R', n, ilo, ihi, work( ileft ),
578 $ work( iright ), n, vsr, ldvsr, ierr )
586 IF( alphai( i ).NE.zero )
THEN
587 IF( ( alphar( i ) / safmax ).GT.( anrmto / anrm ) .OR.
588 $ ( safmin / alphar( i ) ).GT.( anrm / anrmto ) )
THEN
589 work( 1 ) = abs( a( i, i ) / alphar( i ) )
590 beta( i ) = beta( i )*work( 1 )
591 alphar( i ) = alphar( i )*work( 1 )
592 alphai( i ) = alphai( i )*work( 1 )
593 ELSE IF( ( alphai( i ) / safmax ).GT.
594 $ ( anrmto / anrm ) .OR.
595 $ ( safmin / alphai( i ) ).GT.( anrm / anrmto ) )
597 work( 1 ) = abs( a( i, i+1 ) / alphai( i ) )
598 beta( i ) = beta( i )*work( 1 )
599 alphar( i ) = alphar( i )*work( 1 )
600 alphai( i ) = alphai( i )*work( 1 )
608 IF( alphai( i ).NE.zero )
THEN
609 IF( ( beta( i ) / safmax ).GT.( bnrmto / bnrm ) .OR.
610 $ ( safmin / beta( i ) ).GT.( bnrm / bnrmto ) )
THEN
611 work( 1 ) = abs( b( i, i ) / beta( i ) )
612 beta( i ) = beta( i )*work( 1 )
613 alphar( i ) = alphar( i )*work( 1 )
614 alphai( i ) = alphai( i )*work( 1 )
623 CALL
dlascl(
'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
624 CALL
dlascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
625 CALL
dlascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
629 CALL
dlascl(
'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
630 CALL
dlascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
642 cursl = selctg( alphar( i ), alphai( i ), beta( i ) )
643 IF( alphai( i ).EQ.zero )
THEN
647 IF( cursl .AND. .NOT.lastsl )
654 cursl = cursl .OR. lastsl
659 IF( cursl .AND. .NOT.lst2sl )
subroutine dgges(JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK, LWORK, BWORK, INFO)
DGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
subroutine dgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
DGGHRD
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
subroutine dhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
DHGEQZ
subroutine dggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
DGGBAK
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
DGGBAL
subroutine dtgsen(IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO)
DTGSEN