LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ sggsvp3()

subroutine sggsvp3 ( character  JOBU,
character  JOBV,
character  JOBQ,
integer  M,
integer  P,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldb, * )  B,
integer  LDB,
real  TOLA,
real  TOLB,
integer  K,
integer  L,
real, dimension( ldu, * )  U,
integer  LDU,
real, dimension( ldv, * )  V,
integer  LDV,
real, dimension( ldq, * )  Q,
integer  LDQ,
integer, dimension( * )  IWORK,
real, dimension( * )  TAU,
real, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

SGGSVP3

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

Purpose:
 SGGSVP3 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
 SGGSVD3.
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 REAL 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 REAL 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 REAL
[in]TOLB
          TOLB is REAL

          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 REAL 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 REAL 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 REAL 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 REAL array, dimension (N)
[out]WORK
          WORK is REAL 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 SGEQP3 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.

  SGGSVP3 replaces the deprecated subroutine SGGSVP.

Definition at line 269 of file sggsvp3.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  REAL TOLA, TOLB
284 * ..
285 * .. Array Arguments ..
286  INTEGER IWORK( * )
287  REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
288  $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
289 * ..
290 *
291 * =====================================================================
292 *
293 * .. Parameters ..
294  REAL ZERO, ONE
295  parameter( zero = 0.0e+0, one = 1.0e+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 sgeqp3, sgeqr2, sgerq2, slacpy, slapmt,
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 sgeqp3( 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 sgeqp3( m, n, a, lda, iwork, tau, work, -1, info )
366  lwkopt = max( lwkopt, int( work( 1 ) ) )
367  lwkopt = max( 1, lwkopt )
368  work( 1 ) = real( lwkopt )
369  END IF
370 *
371  IF( info.NE.0 ) THEN
372  CALL xerbla( 'SGGSVP3', -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 sgeqp3( p, n, b, ldb, iwork, tau, work, lwork, info )
386 *
387 * Update A := A*P
388 *
389  CALL slapmt( 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 slaset( 'Full', p, p, zero, zero, v, ldv )
404  IF( p.GT.1 )
405  $ CALL slacpy( 'Lower', p-1, n, b( 2, 1 ), ldb, v( 2, 1 ),
406  $ ldv )
407  CALL sorg2r( 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 slaset( '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 slaset( 'Full', n, n, zero, one, q, ldq )
425  CALL slapmt( 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 sgerq2( l, n, b, ldb, tau, work, info )
433 *
434 * Update A := A*Z**T
435 *
436  CALL sormr2( '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 sormr2( 'Right', 'Transpose', n, n, l, b, ldb, tau, q,
444  $ ldq, work, info )
445  END IF
446 *
447 * Clean up B
448 *
449  CALL slaset( '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 sgeqp3( 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 sorm2r( '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 slaset( 'Full', m, m, zero, zero, u, ldu )
489  IF( m.GT.1 )
490  $ CALL slacpy( 'Lower', m-1, n-l, a( 2, 1 ), lda, u( 2, 1 ),
491  $ ldu )
492  CALL sorg2r( 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 slapmt( 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 slaset( '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 sgerq2( 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 sormr2( '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 slaset( '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 sgeqr2( 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 sorm2r( '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 ) = real( lwkopt )
564  RETURN
565 *
566 * End of SGGSVP3
567 *
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine sgerq2(M, N, A, LDA, TAU, WORK, INFO)
SGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm.
Definition: sgerq2.f:123
subroutine sgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
SGEQP3
Definition: sgeqp3.f:151
subroutine sgeqr2(M, N, A, LDA, TAU, WORK, INFO)
SGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
Definition: sgeqr2.f:130
subroutine slapmt(FORWRD, M, N, X, LDX, K)
SLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: slapmt.f:104
subroutine sormr2(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
SORMR2 multiplies a general matrix by the orthogonal matrix from a RQ factorization determined by sge...
Definition: sormr2.f:159
subroutine sorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
SORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition: sorm2r.f:159
subroutine sorg2r(M, N, K, A, LDA, TAU, WORK, INFO)
SORG2R generates all or part of the orthogonal matrix Q from a QR factorization determined by sgeqrf ...
Definition: sorg2r.f:114
Here is the call graph for this function:
Here is the caller graph for this function: