LAPACK 3.12.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 REAL SROUNDUP_LWORK
304 EXTERNAL lsame, sroundup_lwork
305* ..
306* .. External Subroutines ..
307 EXTERNAL sgeqp3, sgeqr2, sgerq2, slacpy, slapmt,
308 $ slaset, sorg2r, sorm2r, sormr2, xerbla
309* ..
310* .. Intrinsic Functions ..
311 INTRINSIC abs, max, min
312* ..
313* .. Executable Statements ..
314*
315* Test the input parameters
316*
317 wantu = lsame( jobu, 'U' )
318 wantv = lsame( jobv, 'V' )
319 wantq = lsame( jobq, 'Q' )
320 forwrd = .true.
321 lquery = ( lwork.EQ.-1 )
322 lwkopt = 1
323*
324* Test the input arguments
325*
326 info = 0
327 IF( .NOT.( wantu .OR. lsame( jobu, 'N' ) ) ) THEN
328 info = -1
329 ELSE IF( .NOT.( wantv .OR. lsame( jobv, 'N' ) ) ) THEN
330 info = -2
331 ELSE IF( .NOT.( wantq .OR. lsame( jobq, 'N' ) ) ) THEN
332 info = -3
333 ELSE IF( m.LT.0 ) THEN
334 info = -4
335 ELSE IF( p.LT.0 ) THEN
336 info = -5
337 ELSE IF( n.LT.0 ) THEN
338 info = -6
339 ELSE IF( lda.LT.max( 1, m ) ) THEN
340 info = -8
341 ELSE IF( ldb.LT.max( 1, p ) ) THEN
342 info = -10
343 ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) ) THEN
344 info = -16
345 ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) ) THEN
346 info = -18
347 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
348 info = -20
349 ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
350 info = -24
351 END IF
352*
353* Compute workspace
354*
355 IF( info.EQ.0 ) THEN
356 CALL sgeqp3( p, n, b, ldb, iwork, tau, work, -1, info )
357 lwkopt = int( work( 1 ) )
358 IF( wantv ) THEN
359 lwkopt = max( lwkopt, p )
360 END IF
361 lwkopt = max( lwkopt, min( n, p ) )
362 lwkopt = max( lwkopt, m )
363 IF( wantq ) THEN
364 lwkopt = max( lwkopt, n )
365 END IF
366 CALL sgeqp3( m, n, a, lda, iwork, tau, work, -1, info )
367 lwkopt = max( lwkopt, int( work( 1 ) ) )
368 lwkopt = max( 1, lwkopt )
369 work( 1 ) = sroundup_lwork( lwkopt )
370 END IF
371*
372 IF( info.NE.0 ) THEN
373 CALL xerbla( 'SGGSVP3', -info )
374 RETURN
375 END IF
376 IF( lquery ) THEN
377 RETURN
378 ENDIF
379*
380* QR with column pivoting of B: B*P = V*( S11 S12 )
381* ( 0 0 )
382*
383 DO 10 i = 1, n
384 iwork( i ) = 0
385 10 CONTINUE
386 CALL sgeqp3( p, n, b, ldb, iwork, tau, work, lwork, info )
387*
388* Update A := A*P
389*
390 CALL slapmt( forwrd, m, n, a, lda, iwork )
391*
392* Determine the effective rank of matrix B.
393*
394 l = 0
395 DO 20 i = 1, min( p, n )
396 IF( abs( b( i, i ) ).GT.tolb )
397 $ l = l + 1
398 20 CONTINUE
399*
400 IF( wantv ) THEN
401*
402* Copy the details of V, and form V.
403*
404 CALL slaset( 'Full', p, p, zero, zero, v, ldv )
405 IF( p.GT.1 )
406 $ CALL slacpy( 'Lower', p-1, n, b( 2, 1 ), ldb, v( 2, 1 ),
407 $ ldv )
408 CALL sorg2r( p, p, min( p, n ), v, ldv, tau, work, info )
409 END IF
410*
411* Clean up B
412*
413 DO 40 j = 1, l - 1
414 DO 30 i = j + 1, l
415 b( i, j ) = zero
416 30 CONTINUE
417 40 CONTINUE
418 IF( p.GT.l )
419 $ CALL slaset( 'Full', p-l, n, zero, zero, b( l+1, 1 ), ldb )
420*
421 IF( wantq ) THEN
422*
423* Set Q = I and Update Q := Q*P
424*
425 CALL slaset( 'Full', n, n, zero, one, q, ldq )
426 CALL slapmt( forwrd, n, n, q, ldq, iwork )
427 END IF
428*
429 IF( p.GE.l .AND. n.NE.l ) THEN
430*
431* RQ factorization of (S11 S12): ( S11 S12 ) = ( 0 S12 )*Z
432*
433 CALL sgerq2( l, n, b, ldb, tau, work, info )
434*
435* Update A := A*Z**T
436*
437 CALL sormr2( 'Right', 'Transpose', m, n, l, b, ldb, tau, a,
438 $ lda, work, info )
439*
440 IF( wantq ) THEN
441*
442* Update Q := Q*Z**T
443*
444 CALL sormr2( 'Right', 'Transpose', n, n, l, b, ldb, tau, q,
445 $ ldq, work, info )
446 END IF
447*
448* Clean up B
449*
450 CALL slaset( 'Full', l, n-l, zero, zero, b, ldb )
451 DO 60 j = n - l + 1, n
452 DO 50 i = j - n + l + 1, l
453 b( i, j ) = zero
454 50 CONTINUE
455 60 CONTINUE
456*
457 END IF
458*
459* Let N-L L
460* A = ( A11 A12 ) M,
461*
462* then the following does the complete QR decomposition of A11:
463*
464* A11 = U*( 0 T12 )*P1**T
465* ( 0 0 )
466*
467 DO 70 i = 1, n - l
468 iwork( i ) = 0
469 70 CONTINUE
470 CALL sgeqp3( m, n-l, a, lda, iwork, tau, work, lwork, info )
471*
472* Determine the effective rank of A11
473*
474 k = 0
475 DO 80 i = 1, min( m, n-l )
476 IF( abs( a( i, i ) ).GT.tola )
477 $ k = k + 1
478 80 CONTINUE
479*
480* Update A12 := U**T*A12, where A12 = A( 1:M, N-L+1:N )
481*
482 CALL sorm2r( 'Left', 'Transpose', m, l, min( m, n-l ), a, lda,
483 $ tau, a( 1, n-l+1 ), lda, work, info )
484*
485 IF( wantu ) THEN
486*
487* Copy the details of U, and form U
488*
489 CALL slaset( 'Full', m, m, zero, zero, u, ldu )
490 IF( m.GT.1 )
491 $ CALL slacpy( 'Lower', m-1, n-l, a( 2, 1 ), lda, u( 2, 1 ),
492 $ ldu )
493 CALL sorg2r( m, m, min( m, n-l ), u, ldu, tau, work, info )
494 END IF
495*
496 IF( wantq ) THEN
497*
498* Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1
499*
500 CALL slapmt( forwrd, n, n-l, q, ldq, iwork )
501 END IF
502*
503* Clean up A: set the strictly lower triangular part of
504* A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0.
505*
506 DO 100 j = 1, k - 1
507 DO 90 i = j + 1, k
508 a( i, j ) = zero
509 90 CONTINUE
510 100 CONTINUE
511 IF( m.GT.k )
512 $ CALL slaset( 'Full', m-k, n-l, zero, zero, a( k+1, 1 ), lda )
513*
514 IF( n-l.GT.k ) THEN
515*
516* RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1
517*
518 CALL sgerq2( k, n-l, a, lda, tau, work, info )
519*
520 IF( wantq ) THEN
521*
522* Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**T
523*
524 CALL sormr2( 'Right', 'Transpose', n, n-l, k, a, lda, tau,
525 $ q, ldq, work, info )
526 END IF
527*
528* Clean up A
529*
530 CALL slaset( 'Full', k, n-l-k, zero, zero, a, lda )
531 DO 120 j = n - l - k + 1, n - l
532 DO 110 i = j - n + l + k + 1, k
533 a( i, j ) = zero
534 110 CONTINUE
535 120 CONTINUE
536*
537 END IF
538*
539 IF( m.GT.k ) THEN
540*
541* QR factorization of A( K+1:M,N-L+1:N )
542*
543 CALL sgeqr2( m-k, l, a( k+1, n-l+1 ), lda, tau, work, info )
544*
545 IF( wantu ) THEN
546*
547* Update U(:,K+1:M) := U(:,K+1:M)*U1
548*
549 CALL sorm2r( 'Right', 'No transpose', m, m-k, min( m-k, l ),
550 $ a( k+1, n-l+1 ), lda, tau, u( 1, k+1 ), ldu,
551 $ work, info )
552 END IF
553*
554* Clean up
555*
556 DO 140 j = n - l + 1, n
557 DO 130 i = j - n + k + l + 1, m
558 a( i, j ) = zero
559 130 CONTINUE
560 140 CONTINUE
561*
562 END IF
563*
564 work( 1 ) = sroundup_lwork( lwkopt )
565 RETURN
566*
567* End of SGGSVP3
568*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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
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
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
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
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
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sroundup_lwork
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
Definition sroundup_lwork.f:59
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
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
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
Here is the call graph for this function:
Here is the caller graph for this function: