LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
August 2015
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 274 of file sggsvp3.f.

274 *
275 * -- LAPACK computational routine (version 3.6.1) --
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  REAL tola, tolb
287 * ..
288 * .. Array Arguments ..
289  INTEGER iwork( * )
290  REAL a( lda, * ), b( ldb, * ), q( ldq, * ),
291  $ tau( * ), u( ldu, * ), v( ldv, * ), work( * )
292 * ..
293 *
294 * =====================================================================
295 *
296 * .. Parameters ..
297  REAL zero, one
298  parameter ( zero = 0.0e+0, one = 1.0e+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 sgeqp3, sgeqr2, sgerq2, slacpy, slapmt,
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 sgeqp3( 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 sgeqp3( m, n, a, lda, iwork, tau, work, -1, info )
369  lwkopt = max( lwkopt, int( work( 1 ) ) )
370  lwkopt = max( 1, lwkopt )
371  work( 1 ) = REAL( lwkopt )
372  END IF
373 *
374  IF( info.NE.0 ) THEN
375  CALL xerbla( 'SGGSVP3', -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 sgeqp3( p, n, b, ldb, iwork, tau, work, lwork, info )
389 *
390 * Update A := A*P
391 *
392  CALL slapmt( 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 slaset( 'Full', p, p, zero, zero, v, ldv )
407  IF( p.GT.1 )
408  $ CALL slacpy( 'Lower', p-1, n, b( 2, 1 ), ldb, v( 2, 1 ),
409  $ ldv )
410  CALL sorg2r( 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 slaset( '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 slaset( 'Full', n, n, zero, one, q, ldq )
428  CALL slapmt( 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 sgerq2( l, n, b, ldb, tau, work, info )
436 *
437 * Update A := A*Z**T
438 *
439  CALL sormr2( '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 sormr2( 'Right', 'Transpose', n, n, l, b, ldb, tau, q,
447  $ ldq, work, info )
448  END IF
449 *
450 * Clean up B
451 *
452  CALL slaset( '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 sgeqp3( 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 sorm2r( '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 slaset( 'Full', m, m, zero, zero, u, ldu )
492  IF( m.GT.1 )
493  $ CALL slacpy( 'Lower', m-1, n-l, a( 2, 1 ), lda, u( 2, 1 ),
494  $ ldu )
495  CALL sorg2r( 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 slapmt( 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 slaset( '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 sgerq2( 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 sormr2( '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 slaset( '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 sgeqr2( 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 sorm2r( '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 ) = REAL( lwkopt )
567  RETURN
568 *
569 * End of SGGSVP3
570 *
subroutine slapmt(FORWRD, M, N, X, LDX, K)
SLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: slapmt.f:106
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
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:123
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:112
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:116
subroutine sgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
SGEQP3
Definition: sgeqp3.f:153
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:161
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:125
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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:161

Here is the call graph for this function:

Here is the caller graph for this function: