LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
August 2015
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 274 of file dggsvp3.f.

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