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 sggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
SGGBAL
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
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 stgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STGEVC
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 sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine shgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SHGEQZ
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 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