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
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
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 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 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 sggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
SGGBAL
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
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 xerbla(SRNAME, INFO)
XERBLA
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
logical function lsame(CA, CB)
LSAME
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.
real function slamch(CMACH)
SLAMCH
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
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 slabad(SMALL, LARGE)
SLABAD
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
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