LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dggsvp3()

subroutine dggsvp3 ( character  jobu,
character  jobv,
character  jobq,
integer  m,
integer  p,
integer  n,
double precision, dimension( lda, * )  a,
integer  lda,
double precision, dimension( ldb, * )  b,
integer  ldb,
double precision  tola,
double precision  tolb,
integer  k,
integer  l,
double precision, dimension( ldu, * )  u,
integer  ldu,
double precision, dimension( ldv, * )  v,
integer  ldv,
double precision, dimension( ldq, * )  q,
integer  ldq,
integer, dimension( * )  iwork,
double precision, dimension( * )  tau,
double precision, dimension( * )  work,
integer  lwork,
integer  info 
)

DGGSVP3

Download DGGSVP3 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DGGSVP3 computes orthogonal matrices U, V and Q such that

                    N-K-L  K    L
  U**T*A*Q =     K ( 0    A12  A13 )  if M-K-L >= 0;
                 L ( 0     0   A23 )
             M-K-L ( 0     0    0  )

                  N-K-L  K    L
         =     K ( 0    A12  A13 )  if M-K-L < 0;
             M-K ( 0     0   A23 )

                  N-K-L  K    L
  V**T*B*Q =   L ( 0     0   B13 )
             P-L ( 0     0    0  )

 where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
 upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
 otherwise A23 is (M-K)-by-L upper trapezoidal.  K+L = the effective
 numerical rank of the (M+P)-by-N matrix (A**T,B**T)**T.

 This decomposition is the preprocessing step for computing the
 Generalized Singular Value Decomposition (GSVD), see subroutine
 DGGSVD3.
Parameters
[in]JOBU
          JOBU is CHARACTER*1
          = 'U':  Orthogonal matrix U is computed;
          = 'N':  U is not computed.
[in]JOBV
          JOBV is CHARACTER*1
          = 'V':  Orthogonal matrix V is computed;
          = 'N':  V is not computed.
[in]JOBQ
          JOBQ is CHARACTER*1
          = 'Q':  Orthogonal matrix Q is computed;
          = 'N':  Q is not computed.
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]P
          P is INTEGER
          The number of rows of the matrix B.  P >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrices A and B.  N >= 0.
[in,out]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, A contains the triangular (or trapezoidal) matrix
          described in the Purpose section.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1,M).
[in,out]B
          B is DOUBLE PRECISION array, dimension (LDB,N)
          On entry, the P-by-N matrix B.
          On exit, B contains the triangular matrix described in
          the Purpose section.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,P).
[in]TOLA
          TOLA is DOUBLE PRECISION
[in]TOLB
          TOLB is DOUBLE PRECISION

          TOLA and TOLB are the thresholds to determine the effective
          numerical rank of matrix B and a subblock of A. Generally,
          they are set to
             TOLA = MAX(M,N)*norm(A)*MACHEPS,
             TOLB = MAX(P,N)*norm(B)*MACHEPS.
          The size of TOLA and TOLB may affect the size of backward
          errors of the decomposition.
[out]K
          K is INTEGER
[out]L
          L is INTEGER

          On exit, K and L specify the dimension of the subblocks
          described in Purpose section.
          K + L = effective numerical rank of (A**T,B**T)**T.
[out]U
          U is DOUBLE PRECISION array, dimension (LDU,M)
          If JOBU = 'U', U contains the orthogonal matrix U.
          If JOBU = 'N', U is not referenced.
[in]LDU
          LDU is INTEGER
          The leading dimension of the array U. LDU >= max(1,M) if
          JOBU = 'U'; LDU >= 1 otherwise.
[out]V
          V is DOUBLE PRECISION array, dimension (LDV,P)
          If JOBV = 'V', V contains the orthogonal matrix V.
          If JOBV = 'N', V is not referenced.
[in]LDV
          LDV is INTEGER
          The leading dimension of the array V. LDV >= max(1,P) if
          JOBV = 'V'; LDV >= 1 otherwise.
[out]Q
          Q is DOUBLE PRECISION array, dimension (LDQ,N)
          If JOBQ = 'Q', Q contains the orthogonal matrix Q.
          If JOBQ = 'N', Q is not referenced.
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q. LDQ >= max(1,N) if
          JOBQ = 'Q'; LDQ >= 1 otherwise.
[out]IWORK
          IWORK is INTEGER array, dimension (N)
[out]TAU
          TAU is DOUBLE PRECISION array, dimension (N)
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  The subroutine uses LAPACK subroutine DGEQP3 for the QR factorization
  with column pivoting to detect the effective numerical rank of the
  a matrix. It may be replaced by a better rank determination strategy.

  DGGSVP3 replaces the deprecated subroutine DGGSVP.

Definition at line 269 of file dggsvp3.f.

272*
273* -- LAPACK computational routine --
274* -- LAPACK is a software package provided by Univ. of Tennessee, --
275* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
276*
277 IMPLICIT NONE
278*
279* .. Scalar Arguments ..
280 CHARACTER JOBQ, JOBU, JOBV
281 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P,
282 $ LWORK
283 DOUBLE PRECISION TOLA, TOLB
284* ..
285* .. Array Arguments ..
286 INTEGER IWORK( * )
287 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
288 $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
289* ..
290*
291* =====================================================================
292*
293* .. Parameters ..
294 DOUBLE PRECISION ZERO, ONE
295 parameter( zero = 0.0d+0, one = 1.0d+0 )
296* ..
297* .. Local Scalars ..
298 LOGICAL FORWRD, WANTQ, WANTU, WANTV, LQUERY
299 INTEGER I, J, LWKOPT
300* ..
301* .. External Functions ..
302 LOGICAL LSAME
303 EXTERNAL lsame
304* ..
305* .. External Subroutines ..
306 EXTERNAL dgeqp3, dgeqr2, dgerq2, dlacpy, dlapmt,
307 $ dlaset, dorg2r, dorm2r, dormr2, xerbla
308* ..
309* .. Intrinsic Functions ..
310 INTRINSIC abs, max, min
311* ..
312* .. Executable Statements ..
313*
314* Test the input parameters
315*
316 wantu = lsame( jobu, 'U' )
317 wantv = lsame( jobv, 'V' )
318 wantq = lsame( jobq, 'Q' )
319 forwrd = .true.
320 lquery = ( lwork.EQ.-1 )
321 lwkopt = 1
322*
323* Test the input arguments
324*
325 info = 0
326 IF( .NOT.( wantu .OR. lsame( jobu, 'N' ) ) ) THEN
327 info = -1
328 ELSE IF( .NOT.( wantv .OR. lsame( jobv, 'N' ) ) ) THEN
329 info = -2
330 ELSE IF( .NOT.( wantq .OR. lsame( jobq, 'N' ) ) ) THEN
331 info = -3
332 ELSE IF( m.LT.0 ) THEN
333 info = -4
334 ELSE IF( p.LT.0 ) THEN
335 info = -5
336 ELSE IF( n.LT.0 ) THEN
337 info = -6
338 ELSE IF( lda.LT.max( 1, m ) ) THEN
339 info = -8
340 ELSE IF( ldb.LT.max( 1, p ) ) THEN
341 info = -10
342 ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) ) THEN
343 info = -16
344 ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) ) THEN
345 info = -18
346 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
347 info = -20
348 ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
349 info = -24
350 END IF
351*
352* Compute workspace
353*
354 IF( info.EQ.0 ) THEN
355 CALL dgeqp3( p, n, b, ldb, iwork, tau, work, -1, info )
356 lwkopt = int( work( 1 ) )
357 IF( wantv ) THEN
358 lwkopt = max( lwkopt, p )
359 END IF
360 lwkopt = max( lwkopt, min( n, p ) )
361 lwkopt = max( lwkopt, m )
362 IF( wantq ) THEN
363 lwkopt = max( lwkopt, n )
364 END IF
365 CALL dgeqp3( m, n, a, lda, iwork, tau, work, -1, info )
366 lwkopt = max( lwkopt, int( work( 1 ) ) )
367 lwkopt = max( 1, lwkopt )
368 work( 1 ) = dble( lwkopt )
369 END IF
370*
371 IF( info.NE.0 ) THEN
372 CALL xerbla( 'DGGSVP3', -info )
373 RETURN
374 END IF
375 IF( lquery ) THEN
376 RETURN
377 ENDIF
378*
379* QR with column pivoting of B: B*P = V*( S11 S12 )
380* ( 0 0 )
381*
382 DO 10 i = 1, n
383 iwork( i ) = 0
384 10 CONTINUE
385 CALL dgeqp3( p, n, b, ldb, iwork, tau, work, lwork, info )
386*
387* Update A := A*P
388*
389 CALL dlapmt( forwrd, m, n, a, lda, iwork )
390*
391* Determine the effective rank of matrix B.
392*
393 l = 0
394 DO 20 i = 1, min( p, n )
395 IF( abs( b( i, i ) ).GT.tolb )
396 $ l = l + 1
397 20 CONTINUE
398*
399 IF( wantv ) THEN
400*
401* Copy the details of V, and form V.
402*
403 CALL dlaset( 'Full', p, p, zero, zero, v, ldv )
404 IF( p.GT.1 )
405 $ CALL dlacpy( 'Lower', p-1, n, b( 2, 1 ), ldb, v( 2, 1 ),
406 $ ldv )
407 CALL dorg2r( p, p, min( p, n ), v, ldv, tau, work, info )
408 END IF
409*
410* Clean up B
411*
412 DO 40 j = 1, l - 1
413 DO 30 i = j + 1, l
414 b( i, j ) = zero
415 30 CONTINUE
416 40 CONTINUE
417 IF( p.GT.l )
418 $ CALL dlaset( 'Full', p-l, n, zero, zero, b( l+1, 1 ), ldb )
419*
420 IF( wantq ) THEN
421*
422* Set Q = I and Update Q := Q*P
423*
424 CALL dlaset( 'Full', n, n, zero, one, q, ldq )
425 CALL dlapmt( forwrd, n, n, q, ldq, iwork )
426 END IF
427*
428 IF( p.GE.l .AND. n.NE.l ) THEN
429*
430* RQ factorization of (S11 S12): ( S11 S12 ) = ( 0 S12 )*Z
431*
432 CALL dgerq2( l, n, b, ldb, tau, work, info )
433*
434* Update A := A*Z**T
435*
436 CALL dormr2( 'Right', 'Transpose', m, n, l, b, ldb, tau, a,
437 $ lda, work, info )
438*
439 IF( wantq ) THEN
440*
441* Update Q := Q*Z**T
442*
443 CALL dormr2( 'Right', 'Transpose', n, n, l, b, ldb, tau, q,
444 $ ldq, work, info )
445 END IF
446*
447* Clean up B
448*
449 CALL dlaset( 'Full', l, n-l, zero, zero, b, ldb )
450 DO 60 j = n - l + 1, n
451 DO 50 i = j - n + l + 1, l
452 b( i, j ) = zero
453 50 CONTINUE
454 60 CONTINUE
455*
456 END IF
457*
458* Let N-L L
459* A = ( A11 A12 ) M,
460*
461* then the following does the complete QR decomposition of A11:
462*
463* A11 = U*( 0 T12 )*P1**T
464* ( 0 0 )
465*
466 DO 70 i = 1, n - l
467 iwork( i ) = 0
468 70 CONTINUE
469 CALL dgeqp3( m, n-l, a, lda, iwork, tau, work, lwork, info )
470*
471* Determine the effective rank of A11
472*
473 k = 0
474 DO 80 i = 1, min( m, n-l )
475 IF( abs( a( i, i ) ).GT.tola )
476 $ k = k + 1
477 80 CONTINUE
478*
479* Update A12 := U**T*A12, where A12 = A( 1:M, N-L+1:N )
480*
481 CALL dorm2r( 'Left', 'Transpose', m, l, min( m, n-l ), a, lda,
482 $ tau, a( 1, n-l+1 ), lda, work, info )
483*
484 IF( wantu ) THEN
485*
486* Copy the details of U, and form U
487*
488 CALL dlaset( 'Full', m, m, zero, zero, u, ldu )
489 IF( m.GT.1 )
490 $ CALL dlacpy( 'Lower', m-1, n-l, a( 2, 1 ), lda, u( 2, 1 ),
491 $ ldu )
492 CALL dorg2r( m, m, min( m, n-l ), u, ldu, tau, work, info )
493 END IF
494*
495 IF( wantq ) THEN
496*
497* Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1
498*
499 CALL dlapmt( forwrd, n, n-l, q, ldq, iwork )
500 END IF
501*
502* Clean up A: set the strictly lower triangular part of
503* A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0.
504*
505 DO 100 j = 1, k - 1
506 DO 90 i = j + 1, k
507 a( i, j ) = zero
508 90 CONTINUE
509 100 CONTINUE
510 IF( m.GT.k )
511 $ CALL dlaset( 'Full', m-k, n-l, zero, zero, a( k+1, 1 ), lda )
512*
513 IF( n-l.GT.k ) THEN
514*
515* RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1
516*
517 CALL dgerq2( k, n-l, a, lda, tau, work, info )
518*
519 IF( wantq ) THEN
520*
521* Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**T
522*
523 CALL dormr2( 'Right', 'Transpose', n, n-l, k, a, lda, tau,
524 $ q, ldq, work, info )
525 END IF
526*
527* Clean up A
528*
529 CALL dlaset( 'Full', k, n-l-k, zero, zero, a, lda )
530 DO 120 j = n - l - k + 1, n - l
531 DO 110 i = j - n + l + k + 1, k
532 a( i, j ) = zero
533 110 CONTINUE
534 120 CONTINUE
535*
536 END IF
537*
538 IF( m.GT.k ) THEN
539*
540* QR factorization of A( K+1:M,N-L+1:N )
541*
542 CALL dgeqr2( m-k, l, a( k+1, n-l+1 ), lda, tau, work, info )
543*
544 IF( wantu ) THEN
545*
546* Update U(:,K+1:M) := U(:,K+1:M)*U1
547*
548 CALL dorm2r( 'Right', 'No transpose', m, m-k, min( m-k, l ),
549 $ a( k+1, n-l+1 ), lda, tau, u( 1, k+1 ), ldu,
550 $ work, info )
551 END IF
552*
553* Clean up
554*
555 DO 140 j = n - l + 1, n
556 DO 130 i = j - n + k + l + 1, m
557 a( i, j ) = zero
558 130 CONTINUE
559 140 CONTINUE
560*
561 END IF
562*
563 work( 1 ) = dble( lwkopt )
564 RETURN
565*
566* End of DGGSVP3
567*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dgeqp3
subroutine dgeqp3(m, n, a, lda, jpvt, tau, work, lwork, info)
DGEQP3
Definition dgeqp3.f:151
dgeqr2
subroutine dgeqr2(m, n, a, lda, tau, work, info)
DGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
Definition dgeqr2.f:130
dgerq2
subroutine dgerq2(m, n, a, lda, tau, work, info)
DGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm.
Definition dgerq2.f:123
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
dlapmt
subroutine dlapmt(forwrd, m, n, x, ldx, k)
DLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition dlapmt.f:104
dlaset
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.
Definition dlaset.f:110
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dorg2r
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 ...
Definition dorg2r.f:114
dorm2r
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...
Definition dorm2r.f:159
dormr2
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...
Definition dormr2.f:159
Here is the call graph for this function:
Here is the caller graph for this function: