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 dlapll(N, X, INCX, Y, INCY, SSMIN)
DLAPLL measures the linear dependence of two vectors.
LOGICAL function lsame(CA, CB)
LSAME
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...
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 dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine xerbla(SRNAME, INFO)
XERBLA
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine dscal(N, DA, DX, INCX)
DSCAL
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
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...