268 SUBROUTINE zgges3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B,
269 $ LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR,
270 $ WORK, LWORK, RWORK, BWORK, INFO )
278 CHARACTER JOBVSL, JOBVSR, SORT
279 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
283 DOUBLE PRECISION RWORK( * )
284 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
285 $ beta( * ), vsl( ldvsl, * ), vsr( ldvsr, * ),
296 DOUBLE PRECISION ZERO, ONE
297 PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
298 COMPLEX*16 CZERO, CONE
299 parameter( czero = ( 0.0d0, 0.0d0 ),
300 $ cone = ( 1.0d0, 0.0d0 ) )
303 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
305 INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
306 $ ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, LWKOPT
307 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
312 DOUBLE PRECISION DIF( 2 )
321 DOUBLE PRECISION DLAMCH, ZLANGE
322 EXTERNAL lsame, dlamch, zlange
331 IF( lsame( jobvsl,
'N' ) )
THEN
334 ELSE IF( lsame( jobvsl,
'V' ) )
THEN
342 IF( lsame( jobvsr,
'N' ) )
THEN
345 ELSE IF( lsame( jobvsr,
'V' ) )
THEN
353 wantst = lsame( sort,
'S' )
358 lquery = ( lwork.EQ.-1 )
359 IF( ijobvl.LE.0 )
THEN
361 ELSE IF( ijobvr.LE.0 )
THEN
363 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort,
'N' ) ) )
THEN
365 ELSE IF( n.LT.0 )
THEN
367 ELSE IF( lda.LT.max( 1, n ) )
THEN
369 ELSE IF( ldb.LT.max( 1, n ) )
THEN
371 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) )
THEN
373 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) )
THEN
375 ELSE IF( lwork.LT.max( 1, 2*n ) .AND. .NOT.lquery )
THEN
382 CALL zgeqrf( n, n, b, ldb, work, work, -1, ierr )
383 lwkopt = max( 1, n + int( work( 1 ) ) )
384 CALL zunmqr(
'L',
'C', n, n, n, b, ldb, work, a, lda, work,
386 lwkopt = max( lwkopt, n + int( work( 1 ) ) )
388 CALL zungqr( n, n, n, vsl, ldvsl, work, work, -1, ierr )
389 lwkopt = max( lwkopt, n + int( work( 1 ) ) )
391 CALL zgghd3( jobvsl, jobvsr, n, 1, n, a, lda, b, ldb, vsl,
392 $ ldvsl, vsr, ldvsr, work, -1, ierr )
393 lwkopt = max( lwkopt, n + int( work( 1 ) ) )
394 CALL zhgeqz(
'S', jobvsl, jobvsr, n, 1, n, a, lda, b, ldb,
395 $ alpha, beta, vsl, ldvsl, vsr, ldvsr, work, -1,
397 lwkopt = max( lwkopt, int( work( 1 ) ) )
399 CALL ztgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
400 $ alpha, beta, vsl, ldvsl, vsr, ldvsr, sdim,
401 $ pvsl, pvsr, dif, work, -1, idum, 1, ierr )
402 lwkopt = max( lwkopt, int( work( 1 ) ) )
404 work( 1 ) = dcmplx( lwkopt )
408 CALL xerbla(
'ZGGES3 ', -info )
410 ELSE IF( lquery )
THEN
424 smlnum = dlamch(
'S' )
425 bignum = one / smlnum
426 CALL dlabad( smlnum, bignum )
427 smlnum = sqrt( smlnum ) / eps
428 bignum = one / smlnum
432 anrm = zlange(
'M', n, n, a, lda, rwork )
434 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
437 ELSE IF( anrm.GT.bignum )
THEN
443 $
CALL zlascl(
'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
447 bnrm = zlange(
'M', n, n, b, ldb, rwork )
449 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
452 ELSE IF( bnrm.GT.bignum )
THEN
458 $
CALL zlascl(
'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
465 CALL zggbal(
'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
466 $ rwork( iright ), rwork( irwrk ), ierr )
470 irows = ihi + 1 - ilo
474 CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
475 $ work( iwrk ), lwork+1-iwrk, ierr )
479 CALL zunmqr(
'L',
'C', irows, icols, irows, b( ilo, ilo ), ldb,
480 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
481 $ lwork+1-iwrk, ierr )
486 CALL zlaset(
'Full', n, n, czero, cone, vsl, ldvsl )
487 IF( irows.GT.1 )
THEN
488 CALL zlacpy(
'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
489 $ vsl( ilo+1, ilo ), ldvsl )
491 CALL zungqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
492 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
498 $
CALL zlaset(
'Full', n, n, czero, cone, vsr, ldvsr )
502 CALL zgghd3( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
503 $ ldvsl, vsr, ldvsr, work( iwrk ), lwork+1-iwrk, ierr )
510 CALL zhgeqz(
'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
511 $ alpha, beta, vsl, ldvsl, vsr, ldvsr, work( iwrk ),
512 $ lwork+1-iwrk, rwork( irwrk ), ierr )
514 IF( ierr.GT.0 .AND. ierr.LE.n )
THEN
516 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n )
THEN
531 $
CALL zlascl(
'G', 0, 0, anrm, anrmto, n, 1, alpha, n, ierr )
533 $
CALL zlascl(
'G', 0, 0, bnrm, bnrmto, n, 1, beta, n, ierr )
538 bwork( i ) = selctg( alpha( i ), beta( i ) )
541 CALL ztgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb, alpha,
542 $ beta, vsl, ldvsl, vsr, ldvsr, sdim, pvsl, pvsr,
543 $ dif, work( iwrk ), lwork-iwrk+1, idum, 1, ierr )
552 $
CALL zggbak(
'P',
'L', n, ilo, ihi, rwork( ileft ),
553 $ rwork( iright ), n, vsl, ldvsl, ierr )
555 $
CALL zggbak(
'P',
'R', n, ilo, ihi, rwork( ileft ),
556 $ rwork( iright ), n, vsr, ldvsr, ierr )
561 CALL zlascl(
'U', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
562 CALL zlascl(
'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
566 CALL zlascl(
'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
567 CALL zlascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
577 cursl = selctg( alpha( i ), beta( i ) )
580 IF( cursl .AND. .NOT.lastsl )
589 work( 1 ) = dcmplx( lwkopt )