377 SUBROUTINE dtgsja( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
378 $ ldb, tola, tolb, alpha, beta, u, ldu, v, ldv,
379 $ q, ldq, work, ncycle, info )
387 CHARACTER jobq, jobu, jobv
388 INTEGER info, k, l, lda, ldb, ldq, ldu, ldv, m, n,
390 DOUBLE PRECISION tola, tolb
393 DOUBLE PRECISION a( lda, * ), alpha( * ), b( ldb, * ),
394 $ beta( * ), q( ldq, * ), u( ldu, * ),
395 $ v( ldv, * ), work( * )
402 parameter( maxit = 40 )
403 DOUBLE PRECISION zero, one
404 parameter( zero = 0.0d+0, one = 1.0d+0 )
408 LOGICAL initq, initu, initv, upper, wantq, wantu, wantv
410 DOUBLE PRECISION a1, a2, a3, b1, b2, b3, csq, csu, csv, error,
411 $ gamma, rwk, snq, snu, snv, ssmin
422 INTRINSIC abs, max, min
428 initu =
lsame( jobu,
'I' )
429 wantu = initu .OR.
lsame( jobu,
'U' )
431 initv =
lsame( jobv,
'I' )
432 wantv = initv .OR.
lsame( jobv,
'V' )
434 initq =
lsame( jobq,
'I' )
435 wantq = initq .OR.
lsame( jobq,
'Q' )
438 IF( .NOT.( initu .OR. wantu .OR.
lsame( jobu,
'N' ) ) )
THEN
440 ELSE IF( .NOT.( initv .OR. wantv .OR.
lsame( jobv,
'N' ) ) )
THEN
442 ELSE IF( .NOT.( initq .OR. wantq .OR.
lsame( jobq,
'N' ) ) )
THEN
444 ELSE IF( m.LT.0 )
THEN
446 ELSE IF( p.LT.0 )
THEN
448 ELSE IF( n.LT.0 )
THEN
450 ELSE IF( lda.LT.max( 1, m ) )
THEN
452 ELSE IF( ldb.LT.max( 1, p ) )
THEN
454 ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) )
THEN
456 ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) )
THEN
458 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) )
THEN
462 CALL
xerbla(
'DTGSJA', -info )
469 $ CALL
dlaset(
'Full', m, m, zero, one, u, ldu )
471 $ CALL
dlaset(
'Full', p, p, zero, one, v, ldv )
473 $ CALL
dlaset(
'Full', n, n, zero, one, q, ldq )
478 DO 40 kcycle = 1, maxit
489 $ a1 = a( k+i, n-l+i )
491 $ a3 = a( k+j, n-l+j )
498 $ a2 = a( k+i, n-l+j )
502 $ a2 = a( k+j, n-l+i )
506 CALL
dlags2( upper, a1, a2, a3, b1, b2, b3, csu, snu,
507 $ csv, snv, csq, snq )
512 $ CALL
drot( l, a( k+j, n-l+1 ), lda, a( k+i, n-l+1 ),
517 CALL
drot( l, b( j, n-l+1 ), ldb, b( i, n-l+1 ), ldb,
523 CALL
drot( min( k+l, m ), a( 1, n-l+j ), 1,
524 $ a( 1, n-l+i ), 1, csq, snq )
526 CALL
drot( l, b( 1, n-l+j ), 1, b( 1, n-l+i ), 1, csq,
531 $ a( k+i, n-l+j ) = zero
535 $ a( k+j, n-l+i ) = zero
541 IF( wantu .AND. k+j.LE.m )
542 $ CALL
drot( m, u( 1, k+j ), 1, u( 1, k+i ), 1, csu,
546 $ CALL
drot( p, v( 1, j ), 1, v( 1, i ), 1, csv, snv )
549 $ CALL
drot( n, q( 1, n-l+j ), 1, q( 1, n-l+i ), 1, csq,
555 IF( .NOT.upper )
THEN
564 DO 30 i = 1, min( l, m-k )
565 CALL
dcopy( l-i+1, a( k+i, n-l+i ), lda, work, 1 )
566 CALL
dcopy( l-i+1, b( i, n-l+i ), ldb, work( l+1 ), 1 )
567 CALL
dlapll( l-i+1, work, 1, work( l+1 ), 1, ssmin )
568 error = max( error, ssmin )
571 IF( abs( error ).LE.min( tola, tolb ) )
595 DO 70 i = 1, min( l, m-k )
600 IF( a1.NE.zero )
THEN
605 IF( gamma.LT.zero )
THEN
606 CALL
dscal( l-i+1, -one, b( i, n-l+i ), ldb )
608 $ CALL
dscal( p, -one, v( 1, i ), 1 )
611 CALL
dlartg( abs( gamma ), one, beta( k+i ), alpha( k+i ),
614 IF( alpha( k+i ).GE.beta( k+i ) )
THEN
615 CALL
dscal( l-i+1, one / alpha( k+i ), a( k+i, n-l+i ),
618 CALL
dscal( l-i+1, one / beta( k+i ), b( i, n-l+i ),
620 CALL
dcopy( l-i+1, b( i, n-l+i ), ldb, a( k+i, n-l+i ),
628 CALL
dcopy( l-i+1, b( i, n-l+i ), ldb, a( k+i, n-l+i ),
637 DO 80 i = m + 1, k + l
643 DO 90 i = k + l + 1, n