271 SUBROUTINE dggsvp3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
272 $ tola, tolb, k, l, u, ldu, v, ldv, q, ldq,
273 $ iwork, tau, work, lwork, info )
283 CHARACTER JOBQ, JOBU, JOBV
284 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P,
286 DOUBLE PRECISION TOLA, TOLB
290 DOUBLE PRECISION A( lda, * ), B( ldb, * ), Q( ldq, * ),
291 $ tau( * ), u( ldu, * ), v( ldv, * ), work( * )
297 DOUBLE PRECISION ZERO, ONE
298 parameter ( zero = 0.0d+0, one = 1.0d+0 )
301 LOGICAL FORWRD, WANTQ, WANTU, WANTV, LQUERY
313 INTRINSIC abs, max, min
319 wantu = lsame( jobu,
'U' )
320 wantv = lsame( jobv,
'V' )
321 wantq = lsame( jobq,
'Q' )
323 lquery = ( lwork.EQ.-1 )
329 IF( .NOT.( wantu .OR. lsame( jobu,
'N' ) ) )
THEN
331 ELSE IF( .NOT.( wantv .OR. lsame( jobv,
'N' ) ) )
THEN
333 ELSE IF( .NOT.( wantq .OR. lsame( jobq,
'N' ) ) )
THEN
335 ELSE IF( m.LT.0 )
THEN
337 ELSE IF( p.LT.0 )
THEN
339 ELSE IF( n.LT.0 )
THEN
341 ELSE IF( lda.LT.max( 1, m ) )
THEN
343 ELSE IF( ldb.LT.max( 1, p ) )
THEN
345 ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) )
THEN
347 ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) )
THEN
349 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) )
THEN
351 ELSE IF( lwork.LT.1 .AND. .NOT.lquery )
THEN
358 CALL dgeqp3( p, n, b, ldb, iwork, tau, work, -1, info )
359 lwkopt = int( work( 1 ) )
361 lwkopt = max( lwkopt, p )
363 lwkopt = max( lwkopt, min( n, p ) )
364 lwkopt = max( lwkopt, m )
366 lwkopt = max( lwkopt, n )
368 CALL dgeqp3( m, n, a, lda, iwork, tau, work, -1, info )
369 lwkopt = max( lwkopt, int( work( 1 ) ) )
370 lwkopt = max( 1, lwkopt )
371 work( 1 ) = dble( lwkopt )
375 CALL xerbla(
'DGGSVP3', -info )
388 CALL dgeqp3( p, n, b, ldb, iwork, tau, work, lwork, info )
392 CALL dlapmt( forwrd, m, n, a, lda, iwork )
397 DO 20 i = 1, min( p, n )
398 IF( abs( b( i, i ) ).GT.tolb )
406 CALL dlaset(
'Full', p, p, zero, zero, v, ldv )
408 $
CALL dlacpy(
'Lower', p-1, n, b( 2, 1 ), ldb, v( 2, 1 ),
410 CALL dorg2r( p, p, min( p, n ), v, ldv, tau, work, info )
421 $
CALL dlaset(
'Full', p-l, n, zero, zero, b( l+1, 1 ), ldb )
427 CALL dlaset(
'Full', n, n, zero, one, q, ldq )
428 CALL dlapmt( forwrd, n, n, q, ldq, iwork )
431 IF( p.GE.l .AND. n.NE.l )
THEN
435 CALL dgerq2( l, n, b, ldb, tau, work, info )
439 CALL dormr2(
'Right',
'Transpose', m, n, l, b, ldb, tau, a,
446 CALL dormr2(
'Right',
'Transpose', n, n, l, b, ldb, tau, q,
452 CALL dlaset(
'Full', l, n-l, zero, zero, b, ldb )
453 DO 60 j = n - l + 1, n
454 DO 50 i = j - n + l + 1, l
472 CALL dgeqp3( m, n-l, a, lda, iwork, tau, work, lwork, info )
477 DO 80 i = 1, min( m, n-l )
478 IF( abs( a( i, i ) ).GT.tola )
484 CALL dorm2r(
'Left',
'Transpose', m, l, min( m, n-l ), a, lda,
485 $ tau, a( 1, n-l+1 ), lda, work, info )
491 CALL dlaset(
'Full', m, m, zero, zero, u, ldu )
493 $
CALL dlacpy(
'Lower', m-1, n-l, a( 2, 1 ), lda, u( 2, 1 ),
495 CALL dorg2r( m, m, min( m, n-l ), u, ldu, tau, work, info )
502 CALL dlapmt( forwrd, n, n-l, q, ldq, iwork )
514 $
CALL dlaset(
'Full', m-k, n-l, zero, zero, a( k+1, 1 ), lda )
520 CALL dgerq2( k, n-l, a, lda, tau, work, info )
526 CALL dormr2(
'Right',
'Transpose', n, n-l, k, a, lda, tau,
527 $ q, ldq, work, info )
532 CALL dlaset(
'Full', k, n-l-k, zero, zero, a, lda )
533 DO 120 j = n - l - k + 1, n - l
534 DO 110 i = j - n + l + k + 1, k
545 CALL dgeqr2( m-k, l, a( k+1, n-l+1 ), lda, tau, work, info )
551 CALL dorm2r(
'Right',
'No transpose', m, m-k, min( m-k, l ),
552 $ a( k+1, n-l+1 ), lda, tau, u( 1, k+1 ), ldu,
558 DO 140 j = n - l + 1, n
559 DO 130 i = j - n + k + l + 1, m
566 work( 1 ) = dble( lwkopt )
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 dgerq2(M, N, A, LDA, TAU, WORK, INFO)
DGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm...
subroutine dgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
DGEQP3
subroutine dgeqr2(M, N, A, LDA, TAU, WORK, INFO)
DGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlapmt(FORWRD, M, N, X, LDX, K)
DLAPMT performs a forward or backward permutation of the columns of a matrix.
subroutine xerbla(SRNAME, INFO)
XERBLA
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 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 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 dggsvp3(JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, TAU, WORK, LWORK, INFO)
DGGSVP3