279 SUBROUTINE sgges3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B,
280 $ LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL,
281 $ VSR, LDVSR, WORK, LWORK, BWORK, INFO )
288 CHARACTER JOBVSL, JOBVSR, SORT
289 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
293 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
294 $ b( ldb, * ), beta( * ), vsl( ldvsl, * ),
295 $ vsr( ldvsr, * ), work( * )
306 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
309 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
310 $ LQUERY, LST2SL, WANTST
311 INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
312 $ ILO, IP, IRIGHT, IROWS, ITAU, IWRK, LWKOPT
313 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
314 $ PVSR, SAFMAX, SAFMIN, SMLNUM
328 EXTERNAL lsame, slamch, slange
331 INTRINSIC abs, max, sqrt
337 IF( lsame( jobvsl,
'N' ) )
THEN
340 ELSE IF( lsame( jobvsl,
'V' ) )
THEN
348 IF( lsame( jobvsr,
'N' ) )
THEN
351 ELSE IF( lsame( jobvsr,
'V' ) )
THEN
359 wantst = lsame( sort,
'S' )
364 lquery = ( lwork.EQ.-1 )
365 IF( ijobvl.LE.0 )
THEN
367 ELSE IF( ijobvr.LE.0 )
THEN
369 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort,
'N' ) ) )
THEN
371 ELSE IF( n.LT.0 )
THEN
373 ELSE IF( lda.LT.max( 1, n ) )
THEN
375 ELSE IF( ldb.LT.max( 1, n ) )
THEN
377 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) )
THEN
379 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) )
THEN
381 ELSE IF( lwork.LT.6*n+16 .AND. .NOT.lquery )
THEN
388 CALL sgeqrf( n, n, b, ldb, work, work, -1, ierr )
389 lwkopt = max( 6*n+16, 3*n+int( work( 1 ) ) )
390 CALL sormqr(
'L',
'T', n, n, n, b, ldb, work, a, lda, work,
392 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
394 CALL sorgqr( n, n, n, vsl, ldvsl, work, work, -1, ierr )
395 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
397 CALL sgghd3( jobvsl, jobvsr, n, 1, n, a, lda, b, ldb, vsl,
398 $ ldvsl, vsr, ldvsr, work, -1, ierr )
399 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
400 CALL slaqz0(
'S', jobvsl, jobvsr, n, 1, n, a, lda, b, ldb,
401 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
402 $ work, -1, 0, ierr )
403 lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
405 CALL stgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
406 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
407 $ sdim, pvsl, pvsr, dif, work, -1, idum, 1,
409 lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
415 CALL xerbla(
'SGGES3 ', -info )
417 ELSE IF( lquery )
THEN
431 safmin = slamch(
'S' )
432 safmax = one / safmin
433 CALL slabad( safmin, safmax )
434 smlnum = sqrt( safmin ) / eps
435 bignum = one / smlnum
439 anrm = slange(
'M', n, n, a, lda, work )
441 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
444 ELSE IF( anrm.GT.bignum )
THEN
449 $
CALL slascl(
'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
453 bnrm = slange(
'M', n, n, b, ldb, work )
455 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
458 ELSE IF( bnrm.GT.bignum )
THEN
463 $
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 )
475 irows = ihi + 1 - ilo
479 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
480 $ work( iwrk ), lwork+1-iwrk, ierr )
484 CALL sormqr(
'L',
'T', irows, icols, irows, b( ilo, ilo ), ldb,
485 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
486 $ lwork+1-iwrk, ierr )
491 CALL slaset(
'Full', n, n, zero, one, vsl, ldvsl )
492 IF( irows.GT.1 )
THEN
493 CALL slacpy(
'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
494 $ vsl( ilo+1, ilo ), ldvsl )
496 CALL sorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
497 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
503 $
CALL slaset(
'Full', n, n, zero, one, vsr, ldvsr )
507 CALL sgghd3( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
508 $ ldvsl, vsr, ldvsr, work( iwrk ), lwork+1-iwrk, ierr )
513 CALL slaqz0(
'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
514 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
515 $ work( iwrk ), lwork+1-iwrk, 0, ierr )
517 IF( ierr.GT.0 .AND. ierr.LE.n )
THEN
519 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n )
THEN
535 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
537 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
541 $
CALL slascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
546 bwork( i ) = selctg( alphar( i ), alphai( i ), beta( i ) )
549 CALL stgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb, alphar,
550 $ alphai, beta, vsl, ldvsl, vsr, ldvsr, sdim, pvsl,
551 $ pvsr, dif, work( iwrk ), lwork-iwrk+1, idum, 1,
561 $
CALL sggbak(
'P',
'L', n, ilo, ihi, work( ileft ),
562 $ work( iright ), n, vsl, ldvsl, ierr )
565 $
CALL sggbak(
'P',
'R', n, ilo, ihi, work( ileft ),
566 $ work( iright ), n, vsr, ldvsr, ierr )
574 IF( alphai( i ).NE.zero )
THEN
575 IF( ( alphar( i )/safmax ).GT.( anrmto/anrm ) .OR.
576 $ ( safmin/alphar( i ) ).GT.( anrm/anrmto ) )
THEN
577 work( 1 ) = abs( a( i, i )/alphar( i ) )
578 beta( i ) = beta( i )*work( 1 )
579 alphar( i ) = alphar( i )*work( 1 )
580 alphai( i ) = alphai( i )*work( 1 )
581 ELSE IF( ( alphai( i )/safmax ).GT.( anrmto/anrm ) .OR.
582 $ ( safmin/alphai( i ) ).GT.( anrm/anrmto ) )
THEN
583 work( 1 ) = abs( a( i, i+1 )/alphai( i ) )
584 beta( i ) = beta( i )*work( 1 )
585 alphar( i ) = alphar( i )*work( 1 )
586 alphai( i ) = alphai( i )*work( 1 )
594 IF( alphai( i ).NE.zero )
THEN
595 IF( ( beta( i )/safmax ).GT.( bnrmto/bnrm ) .OR.
596 $ ( safmin/beta( i ) ).GT.( bnrm/bnrmto ) )
THEN
597 work( 1 ) = abs(b( i, i )/beta( i ))
598 beta( i ) = beta( i )*work( 1 )
599 alphar( i ) = alphar( i )*work( 1 )
600 alphai( i ) = alphai( i )*work( 1 )
609 CALL slascl(
'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
610 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
611 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
615 CALL slascl(
'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
616 CALL slascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
628 cursl = selctg( alphar( i ), alphai( i ), beta( i ) )
629 IF( alphai( i ).EQ.zero )
THEN
633 IF( cursl .AND. .NOT.lastsl )
640 cursl = cursl .OR. lastsl
645 IF( cursl .AND. .NOT.lst2sl )
subroutine slabad(SMALL, LARGE)
SLABAD
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine xerbla(SRNAME, INFO)
XERBLA
recursive subroutine slaqz0(WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, REC, INFO)
SLAQZ0
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 sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
subroutine sgges3(JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK, LWORK, BWORK, INFO)
SGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
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
subroutine sgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SGGHD3