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
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dlapll(N, X, INCX, Y, INCY, SSMIN)
DLAPLL measures the linear dependence of two vectors.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dtgsja(JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, NCYCLE, INFO)
DTGSJA
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dlags2(UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV, SNV, CSQ, SNQ)
DLAGS2 computes 2-by-2 orthogonal matrices U, V, and Q, and applies them to matrices A and B such tha...