377 SUBROUTINE stgsja( 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,
393 REAL A( lda, * ), ALPHA( * ), B( ldb, * ),
394 $ beta( * ), q( ldq, * ), u( ldu, * ),
395 $ v( ldv, * ), work( * )
402 parameter ( maxit = 40 )
404 parameter ( zero = 0.0e+0, one = 1.0e+0 )
408 LOGICAL INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV
410 REAL 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(
'STGSJA', -info )
469 $
CALL slaset(
'Full', m, m, zero, one, u, ldu )
471 $
CALL slaset(
'Full', p, p, zero, one, v, ldv )
473 $
CALL slaset(
'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 slags2( upper, a1, a2, a3, b1, b2, b3, csu, snu,
507 $ csv, snv, csq, snq )
512 $
CALL srot( l, a( k+j, n-l+1 ), lda, a( k+i, n-l+1 ),
517 CALL srot( l, b( j, n-l+1 ), ldb, b( i, n-l+1 ), ldb,
523 CALL srot( min( k+l, m ), a( 1, n-l+j ), 1,
524 $ a( 1, n-l+i ), 1, csq, snq )
526 CALL srot( 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 srot( m, u( 1, k+j ), 1, u( 1, k+i ), 1, csu,
546 $
CALL srot( p, v( 1, j ), 1, v( 1, i ), 1, csv, snv )
549 $
CALL srot( 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 scopy( l-i+1, a( k+i, n-l+i ), lda, work, 1 )
566 CALL scopy( l-i+1, b( i, n-l+i ), ldb, work( l+1 ), 1 )
567 CALL slapll( 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 sscal( l-i+1, -one, b( i, n-l+i ), ldb )
608 $
CALL sscal( p, -one, v( 1, i ), 1 )
611 CALL slartg( abs( gamma ), one, beta( k+i ), alpha( k+i ),
614 IF( alpha( k+i ).GE.beta( k+i ) )
THEN
615 CALL sscal( l-i+1, one / alpha( k+i ), a( k+i, n-l+i ),
618 CALL sscal( l-i+1, one / beta( k+i ), b( i, n-l+i ),
620 CALL scopy( l-i+1, b( i, n-l+i ), ldb, a( k+i, n-l+i ),
628 CALL scopy( 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 slapll(N, X, INCX, Y, INCY, SSMIN)
SLAPLL measures the linear dependence of two vectors.
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine stgsja(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)
STGSJA
subroutine slags2(UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV, SNV, CSQ, SNQ)
SLAGS2 computes 2-by-2 orthogonal matrices U, V, and Q, and applies them to matrices A and B such tha...
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY