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

◆ 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,
307 $ slaset, sorg2r, sorm2r, sormr2, 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 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*
slaset
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
slacpy
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
sgerq2
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
sgeqp3
subroutine sgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
SGEQP3
Definition: sgeqp3.f:151
sgeqr2
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
slapmt
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
sormr2
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
sorm2r
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
sorg2r
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: