378 SUBROUTINE ctgsja( 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,
394 REAL alpha( * ), beta( * )
395 COMPLEX a( lda, * ), b( ldb, * ), q( ldq, * ),
396 $ u( ldu, * ), v( ldv, * ), work( * )
403 parameter( maxit = 40 )
405 parameter( zero = 0.0e+0, one = 1.0e+0 )
407 parameter( czero = ( 0.0e+0, 0.0e+0 ),
408 $ cone = ( 1.0e+0, 0.0e+0 ) )
412 LOGICAL initq, initu, initv, upper, wantq, wantu, wantv
414 REAL a1, a3, b1, b3, csq, csu, csv, error, gamma,
416 COMPLEX a2, b2, snq, snu, snv
427 INTRINSIC abs, conjg, max, min, real
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(
'CTGSJA', -info )
474 $ CALL
claset(
'Full', m, m, czero, cone, u, ldu )
476 $ CALL
claset(
'Full', p, p, czero, cone, v, ldv )
478 $ CALL
claset(
'Full', n, n, czero, cone, q, ldq )
483 DO 40 kcycle = 1, maxit
494 $ a1 =
REAL( A( K+I, N-L+I ) )
496 $ a3 =
REAL( A( K+J, N-L+J ) )
498 b1 =
REAL( B( I, N-L+I ) )
499 b3 =
REAL( B( J, N-L+J ) )
503 $ a2 = a( k+i, n-l+j )
507 $ a2 = a( k+j, n-l+i )
511 CALL
clags2( upper, a1, a2, a3, b1, b2, b3, csu, snu,
512 $ csv, snv, csq, snq )
517 $ CALL
crot( l, a( k+j, n-l+1 ), lda, a( k+i, n-l+1 ),
518 $ lda, csu, conjg( snu ) )
522 CALL
crot( l, b( j, n-l+1 ), ldb, b( i, n-l+1 ), ldb,
523 $ csv, conjg( snv ) )
528 CALL
crot( min( k+l, m ), a( 1, n-l+j ), 1,
529 $ a( 1, n-l+i ), 1, csq, snq )
531 CALL
crot( 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 ) =
REAL( A( K+I, N-L+I ) )
549 $ a( k+j, n-l+j ) =
REAL( A( K+J, N-L+J ) )
550 b( i, n-l+i ) =
REAL( B( I, N-L+I ) )
551 b( j, n-l+j ) =
REAL( B( J, N-L+J ) )
555 IF( wantu .AND. k+j.LE.m )
556 $ CALL
crot( m, u( 1, k+j ), 1, u( 1, k+i ), 1, csu,
560 $ CALL
crot( p, v( 1, j ), 1, v( 1, i ), 1, csv, snv )
563 $ CALL
crot( 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
ccopy( l-i+1, a( k+i, n-l+i ), lda, work, 1 )
580 CALL
ccopy( l-i+1, b( i, n-l+i ), ldb, work( l+1 ), 1 )
581 CALL
clapll( 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 =
REAL( A( K+I, N-L+I ) )
612 b1 =
REAL( B( I, N-L+I ) )
614 IF( a1.NE.zero )
THEN
617 IF( gamma.LT.zero )
THEN
618 CALL
csscal( l-i+1, -one, b( i, n-l+i ), ldb )
620 $ CALL
csscal( p, -one, v( 1, i ), 1 )
623 CALL
slartg( abs( gamma ), one, beta( k+i ), alpha( k+i ),
626 IF( alpha( k+i ).GE.beta( k+i ) )
THEN
627 CALL
csscal( l-i+1, one / alpha( k+i ), a( k+i, n-l+i ),
630 CALL
csscal( l-i+1, one / beta( k+i ), b( i, n-l+i ),
632 CALL
ccopy( l-i+1, b( i, n-l+i ), ldb, a( k+i, n-l+i ),
639 CALL
ccopy( l-i+1, b( i, n-l+i ), ldb, a( k+i, n-l+i ),
646 DO 80 i = m + 1, k + l
652 DO 90 i = k + l + 1, n