LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 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

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.```
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.6.1) --
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,
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 *
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:123
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
subroutine cung2r(M, N, K, A, LDA, TAU, WORK, INFO)
CUNG2R
Definition: cung2r.f:116
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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
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
subroutine cgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
CGEQP3
Definition: cgeqp3.f:161
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
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
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
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: