253 SUBROUTINE dggsvp( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
254 $ tola, tolb, k, l, u, ldu, v, ldv, q, ldq,
255 $ iwork, tau, work, info )
263 CHARACTER jobq, jobu, jobv
264 INTEGER info, k, l, lda, ldb, ldq, ldu, ldv, m, n, p
265 DOUBLE PRECISION tola, tolb
269 DOUBLE PRECISION a( lda, * ), b( ldb, * ), q( ldq, * ),
270 $ tau( * ), u( ldu, * ), v( ldv, * ), work( * )
276 DOUBLE PRECISION zero, one
277 parameter( zero = 0.0d+0, one = 1.0d+0 )
280 LOGICAL forwrd, wantq, wantu, wantv
292 INTRINSIC abs, max, min
298 wantu =
lsame( jobu,
'U' )
299 wantv =
lsame( jobv,
'V' )
300 wantq =
lsame( jobq,
'Q' )
304 IF( .NOT.( wantu .OR.
lsame( jobu,
'N' ) ) )
THEN
306 ELSE IF( .NOT.( wantv .OR.
lsame( jobv,
'N' ) ) )
THEN
308 ELSE IF( .NOT.( wantq .OR.
lsame( jobq,
'N' ) ) )
THEN
310 ELSE IF( m.LT.0 )
THEN
312 ELSE IF( p.LT.0 )
THEN
314 ELSE IF( n.LT.0 )
THEN
316 ELSE IF( lda.LT.max( 1, m ) )
THEN
318 ELSE IF( ldb.LT.max( 1, p ) )
THEN
320 ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) )
THEN
322 ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) )
THEN
324 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) )
THEN
328 CALL
xerbla(
'DGGSVP', -info )
338 CALL
dgeqpf( p, n, b, ldb, iwork, tau, work, info )
342 CALL
dlapmt( forwrd, m, n, a, lda, iwork )
347 DO 20 i = 1, min( p, n )
348 IF( abs( b( i, i ) ).GT.tolb )
356 CALL
dlaset(
'Full', p, p, zero, zero, v, ldv )
358 $ CALL
dlacpy(
'Lower', p-1, n, b( 2, 1 ), ldb, v( 2, 1 ),
360 CALL
dorg2r( p, p, min( p, n ), v, ldv, tau, work, info )
371 $ CALL
dlaset(
'Full', p-l, n, zero, zero, b( l+1, 1 ), ldb )
377 CALL
dlaset(
'Full', n, n, zero, one, q, ldq )
378 CALL
dlapmt( forwrd, n, n, q, ldq, iwork )
381 IF( p.GE.l .AND. n.NE.l )
THEN
385 CALL
dgerq2( l, n, b, ldb, tau, work, info )
389 CALL
dormr2(
'Right',
'Transpose', m, n, l, b, ldb, tau, a,
396 CALL
dormr2(
'Right',
'Transpose', n, n, l, b, ldb, tau, q,
402 CALL
dlaset(
'Full', l, n-l, zero, zero, b, ldb )
403 DO 60 j = n - l + 1, n
404 DO 50 i = j - n + l + 1, l
422 CALL
dgeqpf( m, n-l, a, lda, iwork, tau, work, info )
427 DO 80 i = 1, min( m, n-l )
428 IF( abs( a( i, i ) ).GT.tola )
434 CALL
dorm2r(
'Left',
'Transpose', m, l, min( m, n-l ), a, lda,
435 $ tau, a( 1, n-l+1 ), lda, work, info )
441 CALL
dlaset(
'Full', m, m, zero, zero, u, ldu )
443 $ CALL
dlacpy(
'Lower', m-1, n-l, a( 2, 1 ), lda, u( 2, 1 ),
445 CALL
dorg2r( m, m, min( m, n-l ), u, ldu, tau, work, info )
452 CALL
dlapmt( forwrd, n, n-l, q, ldq, iwork )
464 $ CALL
dlaset(
'Full', m-k, n-l, zero, zero, a( k+1, 1 ), lda )
470 CALL
dgerq2( k, n-l, a, lda, tau, work, info )
476 CALL
dormr2(
'Right',
'Transpose', n, n-l, k, a, lda, tau,
477 $ q, ldq, work, info )
482 CALL
dlaset(
'Full', k, n-l-k, zero, zero, a, lda )
483 DO 120 j = n - l - k + 1, n - l
484 DO 110 i = j - n + l + k + 1, k
495 CALL
dgeqr2( m-k, l, a( k+1, n-l+1 ), lda, tau, work, info )
501 CALL
dorm2r(
'Right',
'No transpose', m, m-k, min( m-k, l ),
502 $ a( k+1, n-l+1 ), lda, tau, u( 1, k+1 ), ldu,
508 DO 140 j = n - l + 1, n
509 DO 130 i = j - n + k + l + 1, m