281 SUBROUTINE sgges( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
282 $ SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR,
283 $ LDVSR, WORK, LWORK, BWORK, INFO )
290 CHARACTER JOBVSL, JOBVSR, SORT
291 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
295 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
296 $ b( ldb, * ), beta( * ), vsl( ldvsl, * ),
297 $ vsr( ldvsr, * ), work( * )
308 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
311 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
312 $ LQUERY, LST2SL, WANTST
313 INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
314 $ ILO, IP, IRIGHT, IROWS, ITAU, IWRK, MAXWRK,
316 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
317 $ PVSR, SAFMAX, SAFMIN, SMLNUM
330 REAL SLAMCH, SLANGE, SROUNDUP_LWORK
331 EXTERNAL lsame, ilaenv, slamch, slange, sroundup_lwork
334 INTRINSIC abs, max, sqrt
340 IF( lsame( jobvsl,
'N' ) )
THEN
343 ELSE IF( lsame( jobvsl,
'V' ) )
THEN
351 IF( lsame( jobvsr,
'N' ) )
THEN
354 ELSE IF( lsame( jobvsr,
'V' ) )
THEN
362 wantst = lsame( sort,
'S' )
367 lquery = ( lwork.EQ.-1 )
368 IF( ijobvl.LE.0 )
THEN
370 ELSE IF( ijobvr.LE.0 )
THEN
372 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort,
'N' ) ) )
THEN
374 ELSE IF( n.LT.0 )
THEN
376 ELSE IF( lda.LT.max( 1, n ) )
THEN
378 ELSE IF( ldb.LT.max( 1, n ) )
THEN
380 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) )
THEN
382 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) )
THEN
395 minwrk = max( 8*n, 6*n + 16 )
396 maxwrk = minwrk - n +
397 $ n*ilaenv( 1,
'SGEQRF',
' ', n, 1, n, 0 )
398 maxwrk = max( maxwrk, minwrk - n +
399 $ n*ilaenv( 1,
'SORMQR',
' ', n, 1, n, -1 ) )
401 maxwrk = max( maxwrk, minwrk - n +
402 $ n*ilaenv( 1,
'SORGQR',
' ', n, 1, n, -1 ) )
408 work( 1 ) = sroundup_lwork(maxwrk)
410 IF( lwork.LT.minwrk .AND. .NOT.lquery )
415 CALL xerbla(
'SGGES ', -info )
417 ELSE IF( lquery )
THEN
431 safmin = slamch(
'S' )
432 safmax = one / safmin
433 smlnum = sqrt( safmin ) / eps
434 bignum = one / smlnum
438 anrm = slange(
'M', n, n, a, lda, work )
440 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
443 ELSE IF( anrm.GT.bignum )
THEN
448 $
CALL slascl(
'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
452 bnrm = slange(
'M', n, n, b, ldb, work )
454 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
457 ELSE IF( bnrm.GT.bignum )
THEN
462 $
CALL slascl(
'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
470 CALL sggbal(
'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
471 $ work( iright ), work( iwrk ), ierr )
476 irows = ihi + 1 - ilo
480 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
481 $ work( iwrk ), lwork+1-iwrk, ierr )
486 CALL sormqr(
'L',
'T', irows, icols, irows, b( ilo, ilo ), ldb,
487 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
488 $ lwork+1-iwrk, ierr )
494 CALL slaset(
'Full', n, n, zero, one, vsl, ldvsl )
495 IF( irows.GT.1 )
THEN
496 CALL slacpy(
'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
497 $ vsl( ilo+1, ilo ), ldvsl )
499 CALL sorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
500 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
506 $
CALL slaset(
'Full', n, n, zero, one, vsr, ldvsr )
511 CALL sgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
512 $ ldvsl, vsr, ldvsr, ierr )
518 CALL shgeqz(
'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
519 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
520 $ work( iwrk ), lwork+1-iwrk, ierr )
522 IF( ierr.GT.0 .AND. ierr.LE.n )
THEN
524 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n )
THEN
541 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
543 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
547 $
CALL slascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
552 bwork( i ) = selctg( alphar( i ), alphai( i ), beta( i ) )
555 CALL stgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb, alphar,
556 $ alphai, beta, vsl, ldvsl, vsr, ldvsr, sdim, pvsl,
557 $ pvsr, dif, work( iwrk ), lwork-iwrk+1, idum, 1,
568 $
CALL sggbak(
'P',
'L', n, ilo, ihi, work( ileft ),
569 $ work( iright ), n, vsl, ldvsl, ierr )
572 $
CALL sggbak(
'P',
'R', n, ilo, ihi, work( ileft ),
573 $ work( iright ), n, vsr, ldvsr, ierr )
581 IF( alphai( i ).NE.zero )
THEN
582 IF( ( alphar( i )/safmax ).GT.( anrmto/anrm ) .OR.
583 $ ( safmin/alphar( i ) ).GT.( anrm/anrmto ) )
THEN
584 work( 1 ) = abs( a( i, i )/alphar( i ) )
585 beta( i ) = beta( i )*work( 1 )
586 alphar( i ) = alphar( i )*work( 1 )
587 alphai( i ) = alphai( i )*work( 1 )
588 ELSE IF( ( alphai( i )/safmax ).GT.( anrmto/anrm ) .OR.
589 $ ( safmin/alphai( i ) ).GT.( anrm/anrmto ) )
THEN
590 work( 1 ) = abs( a( i, i+1 )/alphai( i ) )
591 beta( i ) = beta( i )*work( 1 )
592 alphar( i ) = alphar( i )*work( 1 )
593 alphai( i ) = alphai( i )*work( 1 )
601 IF( alphai( i ).NE.zero )
THEN
602 IF( ( beta( i )/safmax ).GT.( bnrmto/bnrm ) .OR.
603 $ ( safmin/beta( i ) ).GT.( bnrm/bnrmto ) )
THEN
604 work( 1 ) = abs(b( i, i )/beta( i ))
605 beta( i ) = beta( i )*work( 1 )
606 alphar( i ) = alphar( i )*work( 1 )
607 alphai( i ) = alphai( i )*work( 1 )
616 CALL slascl(
'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
617 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
618 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
622 CALL slascl(
'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
623 CALL slascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
635 cursl = selctg( alphar( i ), alphai( i ), beta( i ) )
636 IF( alphai( i ).EQ.zero )
THEN
640 IF( cursl .AND. .NOT.lastsl )
647 cursl = cursl .OR. lastsl
652 IF( cursl .AND. .NOT.lst2sl )
669 work( 1 ) = sroundup_lwork(maxwrk)
subroutine xerbla(srname, info)
subroutine sgeqrf(m, n, a, lda, tau, work, lwork, info)
SGEQRF
subroutine sggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
SGGBAK
subroutine sggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
SGGBAL
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 m...
subroutine sgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
SGGHRD
subroutine shgeqz(job, compq, compz, n, ilo, ihi, h, ldh, t, ldt, alphar, alphai, beta, q, ldq, z, ldz, work, lwork, info)
SHGEQZ
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
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 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 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 sorgqr(m, n, k, a, lda, tau, work, lwork, info)
SORGQR
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR