LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ cggsvp3()

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

CGGSVP3

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

Purpose:
 CGGSVP3 computes unitary matrices U, V and Q such that

                    N-K-L  K    L
  U**H*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**H*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**H,B**H)**H.

 This decomposition is the preprocessing step for computing the
 Generalized Singular Value Decomposition (GSVD), see subroutine
 CGGSVD3.
Parameters
[in]JOBU
          JOBU is CHARACTER*1
          = 'U':  Unitary matrix U is computed;
          = 'N':  U is not computed.
[in]JOBV
          JOBV is CHARACTER*1
          = 'V':  Unitary matrix V is computed;
          = 'N':  V is not computed.
[in]JOBQ
          JOBQ is CHARACTER*1
          = 'Q':  Unitary 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 COMPLEX 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 COMPLEX 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**H,B**H)**H.
[out]U
          U is COMPLEX array, dimension (LDU,M)
          If JOBU = 'U', U contains the unitary 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 COMPLEX array, dimension (LDV,P)
          If JOBV = 'V', V contains the unitary 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 COMPLEX array, dimension (LDQ,N)
          If JOBQ = 'Q', Q contains the unitary 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]RWORK
          RWORK is REAL array, dimension (2*N)
[out]TAU
          TAU is COMPLEX array, dimension (N)
[out]WORK
          WORK is COMPLEX 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 CGEQP3 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.

  CGGSVP3 replaces the deprecated subroutine CGGSVP.

Definition at line 280 of file cggsvp3.f.

280 *
281 * -- LAPACK computational routine (version 3.7.0) --
282 * -- LAPACK is a software package provided by Univ. of Tennessee, --
283 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
284 * August 2015
285 *
286  IMPLICIT NONE
287 *
288 * .. Scalar Arguments ..
289  CHARACTER JOBQ, JOBU, JOBV
290  INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P,
291  $ LWORK
292  REAL TOLA, TOLB
293 * ..
294 * .. Array Arguments ..
295  INTEGER IWORK( * )
296  REAL RWORK( * )
297  COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
298  $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
299 * ..
300 *
301 * =====================================================================
302 *
303 * .. Parameters ..
304  COMPLEX CZERO, CONE
305  parameter( czero = ( 0.0e+0, 0.0e+0 ),
306  $ cone = ( 1.0e+0, 0.0e+0 ) )
307 * ..
308 * .. Local Scalars ..
309  LOGICAL FORWRD, WANTQ, WANTU, WANTV, LQUERY
310  INTEGER I, J, LWKOPT
311 * ..
312 * .. External Functions ..
313  LOGICAL LSAME
314  EXTERNAL lsame
315 * ..
316 * .. External Subroutines ..
317  EXTERNAL cgeqp3, cgeqr2, cgerq2, clacpy, clapmt,
318  $ claset, cung2r, cunm2r, cunmr2, xerbla
319 * ..
320 * .. Intrinsic Functions ..
321  INTRINSIC abs, aimag, max, min, real
322 * ..
323 * .. Executable Statements ..
324 *
325 * Test the input parameters
326 *
327  wantu = lsame( jobu, 'U' )
328  wantv = lsame( jobv, 'V' )
329  wantq = lsame( jobq, 'Q' )
330  forwrd = .true.
331  lquery = ( lwork.EQ.-1 )
332  lwkopt = 1
333 *
334 * Test the input arguments
335 *
336  info = 0
337  IF( .NOT.( wantu .OR. lsame( jobu, 'N' ) ) ) THEN
338  info = -1
339  ELSE IF( .NOT.( wantv .OR. lsame( jobv, 'N' ) ) ) THEN
340  info = -2
341  ELSE IF( .NOT.( wantq .OR. lsame( jobq, 'N' ) ) ) THEN
342  info = -3
343  ELSE IF( m.LT.0 ) THEN
344  info = -4
345  ELSE IF( p.LT.0 ) THEN
346  info = -5
347  ELSE IF( n.LT.0 ) THEN
348  info = -6
349  ELSE IF( lda.LT.max( 1, m ) ) THEN
350  info = -8
351  ELSE IF( ldb.LT.max( 1, p ) ) THEN
352  info = -10
353  ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) ) THEN
354  info = -16
355  ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) ) THEN
356  info = -18
357  ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
358  info = -20
359  ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
360  info = -24
361  END IF
362 *
363 * Compute workspace
364 *
365  IF( info.EQ.0 ) THEN
366  CALL cgeqp3( p, n, b, ldb, iwork, tau, work, -1, rwork, info )
367  lwkopt = int( work( 1 ) )
368  IF( wantv ) THEN
369  lwkopt = max( lwkopt, p )
370  END IF
371  lwkopt = max( lwkopt, min( n, p ) )
372  lwkopt = max( lwkopt, m )
373  IF( wantq ) THEN
374  lwkopt = max( lwkopt, n )
375  END IF
376  CALL cgeqp3( m, n, a, lda, iwork, tau, work, -1, rwork, info )
377  lwkopt = max( lwkopt, int( work( 1 ) ) )
378  lwkopt = max( 1, lwkopt )
379  work( 1 ) = cmplx( lwkopt )
380  END IF
381 *
382  IF( info.NE.0 ) THEN
383  CALL xerbla( 'CGGSVP3', -info )
384  RETURN
385  END IF
386  IF( lquery ) THEN
387  RETURN
388  ENDIF
389 *
390 * QR with column pivoting of B: B*P = V*( S11 S12 )
391 * ( 0 0 )
392 *
393  DO 10 i = 1, n
394  iwork( i ) = 0
395  10 CONTINUE
396  CALL cgeqp3( p, n, b, ldb, iwork, tau, work, lwork, rwork, info )
397 *
398 * Update A := A*P
399 *
400  CALL clapmt( forwrd, m, n, a, lda, iwork )
401 *
402 * Determine the effective rank of matrix B.
403 *
404  l = 0
405  DO 20 i = 1, min( p, n )
406  IF( abs( b( i, i ) ).GT.tolb )
407  $ l = l + 1
408  20 CONTINUE
409 *
410  IF( wantv ) THEN
411 *
412 * Copy the details of V, and form V.
413 *
414  CALL claset( 'Full', p, p, czero, czero, v, ldv )
415  IF( p.GT.1 )
416  $ CALL clacpy( 'Lower', p-1, n, b( 2, 1 ), ldb, v( 2, 1 ),
417  $ ldv )
418  CALL cung2r( p, p, min( p, n ), v, ldv, tau, work, info )
419  END IF
420 *
421 * Clean up B
422 *
423  DO 40 j = 1, l - 1
424  DO 30 i = j + 1, l
425  b( i, j ) = czero
426  30 CONTINUE
427  40 CONTINUE
428  IF( p.GT.l )
429  $ CALL claset( 'Full', p-l, n, czero, czero, b( l+1, 1 ), ldb )
430 *
431  IF( wantq ) THEN
432 *
433 * Set Q = I and Update Q := Q*P
434 *
435  CALL claset( 'Full', n, n, czero, cone, q, ldq )
436  CALL clapmt( forwrd, n, n, q, ldq, iwork )
437  END IF
438 *
439  IF( p.GE.l .AND. n.NE.l ) THEN
440 *
441 * RQ factorization of ( S11 S12 ) = ( 0 S12 )*Z
442 *
443  CALL cgerq2( l, n, b, ldb, tau, work, info )
444 *
445 * Update A := A*Z**H
446 *
447  CALL cunmr2( 'Right', 'Conjugate transpose', m, n, l, b, ldb,
448  $ tau, a, lda, work, info )
449  IF( wantq ) THEN
450 *
451 * Update Q := Q*Z**H
452 *
453  CALL cunmr2( 'Right', 'Conjugate transpose', n, n, l, b,
454  $ ldb, tau, q, ldq, work, info )
455  END IF
456 *
457 * Clean up B
458 *
459  CALL claset( 'Full', l, n-l, czero, czero, b, ldb )
460  DO 60 j = n - l + 1, n
461  DO 50 i = j - n + l + 1, l
462  b( i, j ) = czero
463  50 CONTINUE
464  60 CONTINUE
465 *
466  END IF
467 *
468 * Let N-L L
469 * A = ( A11 A12 ) M,
470 *
471 * then the following does the complete QR decomposition of A11:
472 *
473 * A11 = U*( 0 T12 )*P1**H
474 * ( 0 0 )
475 *
476  DO 70 i = 1, n - l
477  iwork( i ) = 0
478  70 CONTINUE
479  CALL cgeqp3( m, n-l, a, lda, iwork, tau, work, lwork, rwork,
480  $ info )
481 *
482 * Determine the effective rank of A11
483 *
484  k = 0
485  DO 80 i = 1, min( m, n-l )
486  IF( abs( a( i, i ) ).GT.tola )
487  $ k = k + 1
488  80 CONTINUE
489 *
490 * Update A12 := U**H*A12, where A12 = A( 1:M, N-L+1:N )
491 *
492  CALL cunm2r( 'Left', 'Conjugate transpose', m, l, min( m, n-l ),
493  $ a, lda, tau, a( 1, n-l+1 ), lda, work, info )
494 *
495  IF( wantu ) THEN
496 *
497 * Copy the details of U, and form U
498 *
499  CALL claset( 'Full', m, m, czero, czero, u, ldu )
500  IF( m.GT.1 )
501  $ CALL clacpy( 'Lower', m-1, n-l, a( 2, 1 ), lda, u( 2, 1 ),
502  $ ldu )
503  CALL cung2r( m, m, min( m, n-l ), u, ldu, tau, work, info )
504  END IF
505 *
506  IF( wantq ) THEN
507 *
508 * Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1
509 *
510  CALL clapmt( forwrd, n, n-l, q, ldq, iwork )
511  END IF
512 *
513 * Clean up A: set the strictly lower triangular part of
514 * A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0.
515 *
516  DO 100 j = 1, k - 1
517  DO 90 i = j + 1, k
518  a( i, j ) = czero
519  90 CONTINUE
520  100 CONTINUE
521  IF( m.GT.k )
522  $ CALL claset( 'Full', m-k, n-l, czero, czero, a( k+1, 1 ), lda )
523 *
524  IF( n-l.GT.k ) THEN
525 *
526 * RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1
527 *
528  CALL cgerq2( k, n-l, a, lda, tau, work, info )
529 *
530  IF( wantq ) THEN
531 *
532 * Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**H
533 *
534  CALL cunmr2( 'Right', 'Conjugate transpose', n, n-l, k, a,
535  $ lda, tau, q, ldq, work, info )
536  END IF
537 *
538 * Clean up A
539 *
540  CALL claset( 'Full', k, n-l-k, czero, czero, a, lda )
541  DO 120 j = n - l - k + 1, n - l
542  DO 110 i = j - n + l + k + 1, k
543  a( i, j ) = czero
544  110 CONTINUE
545  120 CONTINUE
546 *
547  END IF
548 *
549  IF( m.GT.k ) THEN
550 *
551 * QR factorization of A( K+1:M,N-L+1:N )
552 *
553  CALL cgeqr2( m-k, l, a( k+1, n-l+1 ), lda, tau, work, info )
554 *
555  IF( wantu ) THEN
556 *
557 * Update U(:,K+1:M) := U(:,K+1:M)*U1
558 *
559  CALL cunm2r( 'Right', 'No transpose', m, m-k, min( m-k, l ),
560  $ a( k+1, n-l+1 ), lda, tau, u( 1, k+1 ), ldu,
561  $ work, info )
562  END IF
563 *
564 * Clean up
565 *
566  DO 140 j = n - l + 1, n
567  DO 130 i = j - n + k + l + 1, m
568  a( i, j ) = czero
569  130 CONTINUE
570  140 CONTINUE
571 *
572  END IF
573 *
574  work( 1 ) = cmplx( lwkopt )
575  RETURN
576 *
577 * End of CGGSVP3
578 *
Here is the call graph for this function:
Here is the caller graph for this function:
cgerq2
subroutine cgerq2(M, N, A, LDA, TAU, WORK, INFO)
CGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm.
Definition: cgerq2.f:125
cunmr2
subroutine cunmr2(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
CUNMR2 multiplies a general matrix by the unitary matrix from a RQ factorization determined by cgerqf...
Definition: cunmr2.f:161
cgeqp3
subroutine cgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
CGEQP3
Definition: cgeqp3.f:161
clacpy
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
clapmt
subroutine clapmt(FORWRD, M, N, X, LDX, K)
CLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: clapmt.f:106
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
claset
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:108
cgeqr2
subroutine cgeqr2(M, N, A, LDA, TAU, WORK, INFO)
CGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
Definition: cgeqr2.f:132
cung2r
subroutine cung2r(M, N, K, A, LDA, TAU, WORK, INFO)
CUNG2R
Definition: cung2r.f:116
cunm2r
subroutine cunm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
CUNM2R multiplies a general matrix by the unitary matrix from a QR factorization determined by cgeqrf...
Definition: cunm2r.f:161