283 SUBROUTINE sgges( 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 REAL A( lda, * ), ALPHAI( * ), ALPHAR( * ),
299 $ b( ldb, * ), beta( * ), vsl( ldvsl, * ),
300 $ vsr( ldvsr, * ), work( * )
311 parameter ( zero = 0.0e+0, one = 1.0e+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 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
320 $ pvsr, safmax, safmin, smlnum
335 EXTERNAL lsame, ilaenv, slamch, slange
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,
'SGEQRF',
' ', n, 1, n, 0 )
402 maxwrk = max( maxwrk, minwrk - n +
403 $ n*ilaenv( 1,
'SORMQR',
' ', n, 1, n, -1 ) )
405 maxwrk = max( maxwrk, minwrk - n +
406 $ n*ilaenv( 1,
'SORGQR',
' ', n, 1, n, -1 ) )
414 IF( lwork.LT.minwrk .AND. .NOT.lquery )
419 CALL xerbla(
'SGGES ', -info )
421 ELSE IF( lquery )
THEN
435 safmin = slamch(
'S' )
436 safmax = one / safmin
437 CALL slabad( safmin, safmax )
438 smlnum = sqrt( safmin ) / eps
439 bignum = one / smlnum
443 anrm = slange(
'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 slascl(
'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
457 bnrm = slange(
'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 slascl(
'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
475 CALL sggbal(
'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
476 $ work( iright ), work( iwrk ), ierr )
481 irows = ihi + 1 - ilo
485 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
486 $ work( iwrk ), lwork+1-iwrk, ierr )
491 CALL sormqr(
'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 slaset(
'Full', n, n, zero, one, vsl, ldvsl )
500 IF( irows.GT.1 )
THEN
501 CALL slacpy(
'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
502 $ vsl( ilo+1, ilo ), ldvsl )
504 CALL sorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
505 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
511 $
CALL slaset(
'Full', n, n, zero, one, vsr, ldvsr )
516 CALL sgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
517 $ ldvsl, vsr, ldvsr, ierr )
523 CALL shgeqz(
'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 slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
548 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
552 $
CALL slascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
557 bwork( i ) = selctg( alphar( i ), alphai( i ), beta( i ) )
560 CALL stgsen( 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 sggbak(
'P',
'L', n, ilo, ihi, work( ileft ),
574 $ work( iright ), n, vsl, ldvsl, ierr )
577 $
CALL sggbak(
'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.( anrmto/anrm ) .OR.
594 $ ( safmin/alphai( i ) ).GT.( anrm/anrmto ) )
THEN
595 work( 1 ) = abs( a( i, i+1 )/alphai( i ) )
596 beta( i ) = beta( i )*work( 1 )
597 alphar( i ) = alphar( i )*work( 1 )
598 alphai( i ) = alphai( i )*work( 1 )
606 IF( alphai( i ).NE.zero )
THEN
607 IF( ( beta( i )/safmax ).GT.( bnrmto/bnrm ) .OR.
608 $ ( safmin/beta( i ) ).GT.( bnrm/bnrmto ) )
THEN
609 work( 1 ) = abs(b( i, i )/beta( i ))
610 beta( i ) = beta( i )*work( 1 )
611 alphar( i ) = alphar( i )*work( 1 )
612 alphai( i ) = alphai( i )*work( 1 )
621 CALL slascl(
'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
622 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
623 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
627 CALL slascl(
'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
628 CALL slascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
640 cursl = selctg( alphar( i ), alphai( i ), beta( i ) )
641 IF( alphai( i ).EQ.zero )
THEN
645 IF( cursl .AND. .NOT.lastsl )
652 cursl = cursl .OR. lastsl
657 IF( cursl .AND. .NOT.lst2sl )
subroutine sggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
SGGBAL
subroutine stgsen(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)
STGSEN
subroutine sgges(JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK, LWORK, BWORK, INFO)
SGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
subroutine shgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SHGEQZ
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine sggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
SGGBAK
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine sgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
SGGHRD
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.