255 SUBROUTINE dggsvp( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
256 $ tola, tolb, k, l, u, ldu, v, ldv, q, ldq,
257 $ iwork, tau, work, info )
265 CHARACTER JOBQ, JOBU, JOBV
266 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
267 DOUBLE PRECISION TOLA, TOLB
271 DOUBLE PRECISION A( lda, * ), B( ldb, * ), Q( ldq, * ),
272 $ tau( * ), u( ldu, * ), v( ldv, * ), work( * )
278 DOUBLE PRECISION ZERO, ONE
279 parameter ( zero = 0.0d+0, one = 1.0d+0 )
282 LOGICAL FORWRD, WANTQ, WANTU, WANTV
294 INTRINSIC abs, max, min
300 wantu = lsame( jobu,
'U' )
301 wantv = lsame( jobv,
'V' )
302 wantq = lsame( jobq,
'Q' )
306 IF( .NOT.( wantu .OR. lsame( jobu,
'N' ) ) )
THEN
308 ELSE IF( .NOT.( wantv .OR. lsame( jobv,
'N' ) ) )
THEN
310 ELSE IF( .NOT.( wantq .OR. lsame( jobq,
'N' ) ) )
THEN
312 ELSE IF( m.LT.0 )
THEN
314 ELSE IF( p.LT.0 )
THEN
316 ELSE IF( n.LT.0 )
THEN
318 ELSE IF( lda.LT.max( 1, m ) )
THEN
320 ELSE IF( ldb.LT.max( 1, p ) )
THEN
322 ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) )
THEN
324 ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) )
THEN
326 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) )
THEN
330 CALL xerbla(
'DGGSVP', -info )
340 CALL dgeqpf( p, n, b, ldb, iwork, tau, work, info )
344 CALL dlapmt( forwrd, m, n, a, lda, iwork )
349 DO 20 i = 1, min( p, n )
350 IF( abs( b( i, i ) ).GT.tolb )
358 CALL dlaset(
'Full', p, p, zero, zero, v, ldv )
360 $
CALL dlacpy(
'Lower', p-1, n, b( 2, 1 ), ldb, v( 2, 1 ),
362 CALL dorg2r( p, p, min( p, n ), v, ldv, tau, work, info )
373 $
CALL dlaset(
'Full', p-l, n, zero, zero, b( l+1, 1 ), ldb )
379 CALL dlaset(
'Full', n, n, zero, one, q, ldq )
380 CALL dlapmt( forwrd, n, n, q, ldq, iwork )
383 IF( p.GE.l .AND. n.NE.l )
THEN
387 CALL dgerq2( l, n, b, ldb, tau, work, info )
391 CALL dormr2(
'Right',
'Transpose', m, n, l, b, ldb, tau, a,
398 CALL dormr2(
'Right',
'Transpose', n, n, l, b, ldb, tau, q,
404 CALL dlaset(
'Full', l, n-l, zero, zero, b, ldb )
405 DO 60 j = n - l + 1, n
406 DO 50 i = j - n + l + 1, l
424 CALL dgeqpf( m, n-l, a, lda, iwork, tau, work, info )
429 DO 80 i = 1, min( m, n-l )
430 IF( abs( a( i, i ) ).GT.tola )
436 CALL dorm2r(
'Left',
'Transpose', m, l, min( m, n-l ), a, lda,
437 $ tau, a( 1, n-l+1 ), lda, work, info )
443 CALL dlaset(
'Full', m, m, zero, zero, u, ldu )
445 $
CALL dlacpy(
'Lower', m-1, n-l, a( 2, 1 ), lda, u( 2, 1 ),
447 CALL dorg2r( m, m, min( m, n-l ), u, ldu, tau, work, info )
454 CALL dlapmt( forwrd, n, n-l, q, ldq, iwork )
466 $
CALL dlaset(
'Full', m-k, n-l, zero, zero, a( k+1, 1 ), lda )
472 CALL dgerq2( k, n-l, a, lda, tau, work, info )
478 CALL dormr2(
'Right',
'Transpose', n, n-l, k, a, lda, tau,
479 $ q, ldq, work, info )
484 CALL dlaset(
'Full', k, n-l-k, zero, zero, a, lda )
485 DO 120 j = n - l - k + 1, n - l
486 DO 110 i = j - n + l + k + 1, k
497 CALL dgeqr2( m-k, l, a( k+1, n-l+1 ), lda, tau, work, info )
503 CALL dorm2r(
'Right',
'No transpose', m, m-k, min( m-k, l ),
504 $ a( k+1, n-l+1 ), lda, tau, u( 1, k+1 ), ldu,
510 DO 140 j = n - l + 1, n
511 DO 130 i = j - n + k + l + 1, m
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dormr2(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORMR2 multiplies a general matrix by the orthogonal matrix from a RQ factorization determined by sge...
subroutine dgerq2(M, N, A, LDA, TAU, WORK, INFO)
DGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm...
subroutine dgeqr2(M, N, A, LDA, TAU, WORK, INFO)
DGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
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 dorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
subroutine dggsvp(JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, TAU, WORK, INFO)
DGGSVP
subroutine dorg2r(M, N, K, A, LDA, TAU, WORK, INFO)
DORG2R generates all or part of the orthogonal matrix Q from a QR factorization determined by sgeqrf ...
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dgeqpf(M, N, A, LDA, JPVT, TAU, WORK, INFO)
DGEQPF
subroutine dlapmt(FORWRD, M, N, X, LDX, K)
DLAPMT performs a forward or backward permutation of the columns of a matrix.