378 SUBROUTINE ztgsja( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
379 $ LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV,
380 $ Q, LDQ, WORK, NCYCLE, INFO )
388 CHARACTER JOBQ, JOBU, JOBV
389 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N,
391 DOUBLE PRECISION TOLA, TOLB
394 DOUBLE PRECISION ALPHA( * ), BETA( * )
395 COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
396 $ u( ldu, * ), v( ldv, * ), work( * )
403 PARAMETER ( MAXIT = 40 )
404 DOUBLE PRECISION ZERO, ONE
405 parameter( zero = 0.0d+0, one = 1.0d+0 )
406 COMPLEX*16 CZERO, CONE
407 parameter( czero = ( 0.0d+0, 0.0d+0 ),
408 $ cone = ( 1.0d+0, 0.0d+0 ) )
412 LOGICAL INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV
414 DOUBLE PRECISION A1, A3, B1, B3, CSQ, CSU, CSV, ERROR, GAMMA,
416 COMPLEX*16 A2, B2, SNQ, SNU, SNV
427 INTRINSIC abs, dble, dconjg, max, min
433 initu = lsame( jobu,
'I' )
434 wantu = initu .OR. lsame( jobu,
'U' )
436 initv = lsame( jobv,
'I' )
437 wantv = initv .OR. lsame( jobv,
'V' )
439 initq = lsame( jobq,
'I' )
440 wantq = initq .OR. lsame( jobq,
'Q' )
443 IF( .NOT.( initu .OR. wantu .OR. lsame( jobu,
'N' ) ) )
THEN
445 ELSE IF( .NOT.( initv .OR. wantv .OR. lsame( jobv,
'N' ) ) )
THEN
447 ELSE IF( .NOT.( initq .OR. wantq .OR. lsame( jobq,
'N' ) ) )
THEN
449 ELSE IF( m.LT.0 )
THEN
451 ELSE IF( p.LT.0 )
THEN
453 ELSE IF( n.LT.0 )
THEN
455 ELSE IF( lda.LT.max( 1, m ) )
THEN
457 ELSE IF( ldb.LT.max( 1, p ) )
THEN
459 ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) )
THEN
461 ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) )
THEN
463 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) )
THEN
467 CALL xerbla(
'ZTGSJA', -info )
474 $
CALL zlaset(
'Full', m, m, czero, cone, u, ldu )
476 $
CALL zlaset(
'Full', p, p, czero, cone, v, ldv )
478 $
CALL zlaset(
'Full', n, n, czero, cone, q, ldq )
483 DO 40 kcycle = 1, maxit
494 $ a1 = dble( a( k+i, n-l+i ) )
496 $ a3 = dble( a( k+j, n-l+j ) )
498 b1 = dble( b( i, n-l+i ) )
499 b3 = dble( b( j, n-l+j ) )
503 $ a2 = a( k+i, n-l+j )
507 $ a2 = a( k+j, n-l+i )
511 CALL zlags2( upper, a1, a2, a3, b1, b2, b3, csu, snu,
512 $ csv, snv, csq, snq )
517 $
CALL zrot( l, a( k+j, n-l+1 ), lda, a( k+i, n-l+1 ),
518 $ lda, csu, dconjg( snu ) )
522 CALL zrot( l, b( j, n-l+1 ), ldb, b( i, n-l+1 ), ldb,
523 $ csv, dconjg( snv ) )
528 CALL zrot( min( k+l, m ), a( 1, n-l+j ), 1,
529 $ a( 1, n-l+i ), 1, csq, snq )
531 CALL zrot( l, b( 1, n-l+j ), 1, b( 1, n-l+i ), 1, csq,
536 $ a( k+i, n-l+j ) = czero
537 b( i, n-l+j ) = czero
540 $ a( k+j, n-l+i ) = czero
541 b( j, n-l+i ) = czero
547 $ a( k+i, n-l+i ) = dble( a( k+i, n-l+i ) )
549 $ a( k+j, n-l+j ) = dble( a( k+j, n-l+j ) )
550 b( i, n-l+i ) = dble( b( i, n-l+i ) )
551 b( j, n-l+j ) = dble( b( j, n-l+j ) )
555 IF( wantu .AND. k+j.LE.m )
556 $
CALL zrot( m, u( 1, k+j ), 1, u( 1, k+i ), 1, csu,
560 $
CALL zrot( p, v( 1, j ), 1, v( 1, i ), 1, csv, snv )
563 $
CALL zrot( n, q( 1, n-l+j ), 1, q( 1, n-l+i ), 1, csq,
569 IF( .NOT.upper )
THEN
578 DO 30 i = 1, min( l, m-k )
579 CALL zcopy( l-i+1, a( k+i, n-l+i ), lda, work, 1 )
580 CALL zcopy( l-i+1, b( i, n-l+i ), ldb, work( l+1 ), 1 )
581 CALL zlapll( l-i+1, work, 1, work( l+1 ), 1, ssmin )
582 error = max( error, ssmin )
585 IF( abs( error ).LE.min( tola, tolb ) )
609 DO 70 i = 1, min( l, m-k )
611 a1 = dble( a( k+i, n-l+i ) )
612 b1 = dble( b( i, n-l+i ) )
614 IF( a1.NE.zero )
THEN
617 IF( gamma.LT.zero )
THEN
618 CALL zdscal( l-i+1, -one, b( i, n-l+i ), ldb )
620 $
CALL zdscal( p, -one, v( 1, i ), 1 )
623 CALL dlartg( abs( gamma ), one, beta( k+i ), alpha( k+i ),
626 IF( alpha( k+i ).GE.beta( k+i ) )
THEN
627 CALL zdscal( l-i+1, one / alpha( k+i ), a( k+i, n-l+i ),
630 CALL zdscal( l-i+1, one / beta( k+i ), b( i, n-l+i ),
632 CALL zcopy( l-i+1, b( i, n-l+i ), ldb, a( k+i, n-l+i ),
640 CALL zcopy( l-i+1, b( i, n-l+i ), ldb, a( k+i, n-l+i ),
647 DO 80 i = m + 1, k + l
653 DO 90 i = k + l + 1, n