306 SUBROUTINE sgegv( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
307 $ beta, vl, ldvl, vr, ldvr, work, lwork, info )
315 CHARACTER JOBVL, JOBVR
316 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
319 REAL A( lda, * ), ALPHAI( * ), ALPHAR( * ),
320 $ b( ldb, * ), beta( * ), vl( ldvl, * ),
321 $ vr( ldvr, * ), work( * )
328 parameter( zero = 0.0e0, one = 1.0e0 )
331 LOGICAL ILIMIT, ILV, ILVL, ILVR, LQUERY
333 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
334 $ in, iright, irows, itau, iwork, jc, jr, lopt,
335 $ lwkmin, lwkopt, nb, nb1, nb2, nb3
336 REAL ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
337 $ bnrm1, bnrm2, eps, onepls, safmax, safmin,
338 $ salfai, salfar, sbeta, scale, temp
351 EXTERNAL ilaenv, lsame, slamch, slange
354 INTRINSIC abs, int, max
360 IF( lsame( jobvl,
'N' ) )
THEN
363 ELSE IF( lsame( jobvl,
'V' ) )
THEN
371 IF( lsame( jobvr,
'N' ) )
THEN
374 ELSE IF( lsame( jobvr,
'V' ) )
THEN
385 lwkmin = max( 8*n, 1 )
388 lquery = ( lwork.EQ.-1 )
390 IF( ijobvl.LE.0 )
THEN
392 ELSE IF( ijobvr.LE.0 )
THEN
394 ELSE IF( n.LT.0 )
THEN
396 ELSE IF( lda.LT.max( 1, n ) )
THEN
398 ELSE IF( ldb.LT.max( 1, n ) )
THEN
400 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) )
THEN
402 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) )
THEN
404 ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery )
THEN
409 nb1 = ilaenv( 1,
'SGEQRF',
' ', n, n, -1, -1 )
410 nb2 = ilaenv( 1,
'SORMQR',
' ', n, n, n, -1 )
411 nb3 = ilaenv( 1,
'SORGQR',
' ', n, n, n, -1 )
412 nb = max( nb1, nb2, nb3 )
413 lopt = 2*n + max( 6*n, n*(nb+1) )
418 CALL
xerbla(
'SGEGV ', -info )
420 ELSE IF( lquery )
THEN
431 eps = slamch(
'E' )*slamch(
'B' )
432 safmin = slamch(
'S' )
433 safmin = safmin + safmin
434 safmax = one / safmin
435 onepls = one + ( 4*eps )
439 anrm = slange(
'M', n, n, a, lda, work )
442 IF( anrm.LT.one )
THEN
443 IF( safmax*anrm.LT.one )
THEN
449 IF( anrm.GT.zero )
THEN
450 CALL
slascl(
'G', -1, -1, anrm, one, n, n, a, lda, iinfo )
451 IF( iinfo.NE.0 )
THEN
459 bnrm = slange(
'M', n, n, b, ldb, work )
462 IF( bnrm.LT.one )
THEN
463 IF( safmax*bnrm.LT.one )
THEN
469 IF( bnrm.GT.zero )
THEN
470 CALL
slascl(
'G', -1, -1, bnrm, one, n, n, b, ldb, iinfo )
471 IF( iinfo.NE.0 )
THEN
484 CALL
sggbal(
'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
485 $ work( iright ), work( iwork ), iinfo )
486 IF( iinfo.NE.0 )
THEN
495 irows = ihi + 1 - ilo
503 CALL
sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
504 $ work( iwork ), lwork+1-iwork, iinfo )
506 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
507 IF( iinfo.NE.0 )
THEN
512 CALL
sormqr(
'L',
'T', irows, icols, irows, b( ilo, ilo ), ldb,
513 $ work( itau ), a( ilo, ilo ), lda, work( iwork ),
514 $ lwork+1-iwork, iinfo )
516 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
517 IF( iinfo.NE.0 )
THEN
523 CALL
slaset(
'Full', n, n, zero, one, vl, ldvl )
524 CALL
slacpy(
'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
525 $ vl( ilo+1, ilo ), ldvl )
526 CALL
sorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
527 $ work( itau ), work( iwork ), lwork+1-iwork,
530 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
531 IF( iinfo.NE.0 )
THEN
538 $ CALL
slaset(
'Full', n, n, zero, one, vr, ldvr )
546 CALL
sgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
547 $ ldvl, vr, ldvr, iinfo )
549 CALL
sgghrd(
'N',
'N', irows, 1, irows, a( ilo, ilo ), lda,
550 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, iinfo )
552 IF( iinfo.NE.0 )
THEN
567 CALL
shgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
568 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
569 $ work( iwork ), lwork+1-iwork, iinfo )
571 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
572 IF( iinfo.NE.0 )
THEN
573 IF( iinfo.GT.0 .AND. iinfo.LE.n )
THEN
575 ELSE IF( iinfo.GT.n .AND. iinfo.LE.2*n )
THEN
597 CALL
stgevc( chtemp,
'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
598 $ vr, ldvr, n, in, work( iwork ), iinfo )
599 IF( iinfo.NE.0 )
THEN
607 CALL
sggbak(
'P',
'L', n, ilo, ihi, work( ileft ),
608 $ work( iright ), n, vl, ldvl, iinfo )
609 IF( iinfo.NE.0 )
THEN
614 IF( alphai( jc ).LT.zero )
617 IF( alphai( jc ).EQ.zero )
THEN
619 temp = max( temp, abs( vl( jr, jc ) ) )
623 temp = max( temp, abs( vl( jr, jc ) )+
624 $ abs( vl( jr, jc+1 ) ) )
630 IF( alphai( jc ).EQ.zero )
THEN
632 vl( jr, jc ) = vl( jr, jc )*temp
636 vl( jr, jc ) = vl( jr, jc )*temp
637 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
643 CALL
sggbak(
'P',
'R', n, ilo, ihi, work( ileft ),
644 $ work( iright ), n, vr, ldvr, iinfo )
645 IF( iinfo.NE.0 )
THEN
650 IF( alphai( jc ).LT.zero )
653 IF( alphai( jc ).EQ.zero )
THEN
655 temp = max( temp, abs( vr( jr, jc ) ) )
659 temp = max( temp, abs( vr( jr, jc ) )+
660 $ abs( vr( jr, jc+1 ) ) )
666 IF( alphai( jc ).EQ.zero )
THEN
668 vr( jr, jc ) = vr( jr, jc )*temp
672 vr( jr, jc ) = vr( jr, jc )*temp
673 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
692 absar = abs( alphar( jc ) )
693 absai = abs( alphai( jc ) )
694 absb = abs( beta( jc ) )
695 salfar = anrm*alphar( jc )
696 salfai = anrm*alphai( jc )
697 sbeta = bnrm*beta( jc )
703 IF( abs( salfai ).LT.safmin .AND. absai.GE.
704 $ max( safmin, eps*absar, eps*absb ) )
THEN
706 scale = ( onepls*safmin / anrm1 ) /
707 $ max( onepls*safmin, anrm2*absai )
709 ELSE IF( salfai.EQ.zero )
THEN
714 IF( alphai( jc ).LT.zero .AND. jc.GT.1 )
THEN
715 alphai( jc-1 ) = zero
716 ELSE IF( alphai( jc ).GT.zero .AND. jc.LT.n )
THEN
717 alphai( jc+1 ) = zero
723 IF( abs( salfar ).LT.safmin .AND. absar.GE.
724 $ max( safmin, eps*absai, eps*absb ) )
THEN
726 scale = max( scale, ( onepls*safmin / anrm1 ) /
727 $ max( onepls*safmin, anrm2*absar ) )
732 IF( abs( sbeta ).LT.safmin .AND. absb.GE.
733 $ max( safmin, eps*absar, eps*absai ) )
THEN
735 scale = max( scale, ( onepls*safmin / bnrm1 ) /
736 $ max( onepls*safmin, bnrm2*absb ) )
742 temp = ( scale*safmin )*max( abs( salfar ), abs( salfai ),
745 $ scale = scale / temp
753 salfar = ( scale*alphar( jc ) )*anrm
754 salfai = ( scale*alphai( jc ) )*anrm
755 sbeta = ( scale*beta( jc ) )*bnrm
757 alphar( jc ) = salfar
758 alphai( jc ) = salfai
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 sggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
SGGBAL
subroutine sgegv(JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO)
SGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
subroutine sgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
SGGHRD
subroutine stgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STGEVC
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine shgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SHGEQZ
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
subroutine sggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
SGGBAK