LAPACK  3.10.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.
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,
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 *
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
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
DGEQP3
Definition: dgeqp3.f:151
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
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
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
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
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
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: