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

## ◆ 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

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.```
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 275 of file cggsvp3.f.

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