LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
cggsvp.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine cggsvp (JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK, TAU, WORK, INFO)
 CGGSVP More...
 

Function/Subroutine Documentation

subroutine cggsvp ( 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  INFO 
)

CGGSVP

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

Purpose:
 CGGSVP 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
 CGGSVD.
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(3*N,M,P))
[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
November 2011
Further Details:
The subroutine uses LAPACK subroutine CGEQPF 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.

Definition at line 262 of file cggsvp.f.

262 *
263 * -- LAPACK computational routine (version 3.4.0) --
264 * -- LAPACK is a software package provided by Univ. of Tennessee, --
265 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
266 * November 2011
267 *
268 * .. Scalar Arguments ..
269  CHARACTER jobq, jobu, jobv
270  INTEGER info, k, l, lda, ldb, ldq, ldu, ldv, m, n, p
271  REAL tola, tolb
272 * ..
273 * .. Array Arguments ..
274  INTEGER iwork( * )
275  REAL rwork( * )
276  COMPLEX a( lda, * ), b( ldb, * ), q( ldq, * ),
277  $ tau( * ), u( ldu, * ), v( ldv, * ), work( * )
278 * ..
279 *
280 * =====================================================================
281 *
282 * .. Parameters ..
283  COMPLEX czero, cone
284  parameter( czero = ( 0.0e+0, 0.0e+0 ),
285  $ cone = ( 1.0e+0, 0.0e+0 ) )
286 * ..
287 * .. Local Scalars ..
288  LOGICAL forwrd, wantq, wantu, wantv
289  INTEGER i, j
290  COMPLEX t
291 * ..
292 * .. External Functions ..
293  LOGICAL lsame
294  EXTERNAL lsame
295 * ..
296 * .. External Subroutines ..
297  EXTERNAL cgeqpf, cgeqr2, cgerq2, clacpy, clapmt, claset,
299 * ..
300 * .. Intrinsic Functions ..
301  INTRINSIC abs, aimag, max, min, real
302 * ..
303 * .. Statement Functions ..
304  REAL cabs1
305 * ..
306 * .. Statement Function definitions ..
307  cabs1( t ) = abs( REAL( T ) ) + abs( aimag( t ) )
308 * ..
309 * .. Executable Statements ..
310 *
311 * Test the input parameters
312 *
313  wantu = lsame( jobu, 'U' )
314  wantv = lsame( jobv, 'V' )
315  wantq = lsame( jobq, 'Q' )
316  forwrd = .true.
317 *
318  info = 0
319  IF( .NOT.( wantu .OR. lsame( jobu, 'N' ) ) ) THEN
320  info = -1
321  ELSE IF( .NOT.( wantv .OR. lsame( jobv, 'N' ) ) ) THEN
322  info = -2
323  ELSE IF( .NOT.( wantq .OR. lsame( jobq, 'N' ) ) ) THEN
324  info = -3
325  ELSE IF( m.LT.0 ) THEN
326  info = -4
327  ELSE IF( p.LT.0 ) THEN
328  info = -5
329  ELSE IF( n.LT.0 ) THEN
330  info = -6
331  ELSE IF( lda.LT.max( 1, m ) ) THEN
332  info = -8
333  ELSE IF( ldb.LT.max( 1, p ) ) THEN
334  info = -10
335  ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) ) THEN
336  info = -16
337  ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) ) THEN
338  info = -18
339  ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
340  info = -20
341  END IF
342  IF( info.NE.0 ) THEN
343  CALL xerbla( 'CGGSVP', -info )
344  RETURN
345  END IF
346 *
347 * QR with column pivoting of B: B*P = V*( S11 S12 )
348 * ( 0 0 )
349 *
350  DO 10 i = 1, n
351  iwork( i ) = 0
352  10 CONTINUE
353  CALL cgeqpf( p, n, b, ldb, iwork, tau, work, rwork, info )
354 *
355 * Update A := A*P
356 *
357  CALL clapmt( forwrd, m, n, a, lda, iwork )
358 *
359 * Determine the effective rank of matrix B.
360 *
361  l = 0
362  DO 20 i = 1, min( p, n )
363  IF( cabs1( b( i, i ) ).GT.tolb )
364  $ l = l + 1
365  20 CONTINUE
366 *
367  IF( wantv ) THEN
368 *
369 * Copy the details of V, and form V.
370 *
371  CALL claset( 'Full', p, p, czero, czero, v, ldv )
372  IF( p.GT.1 )
373  $ CALL clacpy( 'Lower', p-1, n, b( 2, 1 ), ldb, v( 2, 1 ),
374  $ ldv )
375  CALL cung2r( p, p, min( p, n ), v, ldv, tau, work, info )
376  END IF
377 *
378 * Clean up B
379 *
380  DO 40 j = 1, l - 1
381  DO 30 i = j + 1, l
382  b( i, j ) = czero
383  30 CONTINUE
384  40 CONTINUE
385  IF( p.GT.l )
386  $ CALL claset( 'Full', p-l, n, czero, czero, b( l+1, 1 ), ldb )
387 *
388  IF( wantq ) THEN
389 *
390 * Set Q = I and Update Q := Q*P
391 *
392  CALL claset( 'Full', n, n, czero, cone, q, ldq )
393  CALL clapmt( forwrd, n, n, q, ldq, iwork )
394  END IF
395 *
396  IF( p.GE.l .AND. n.NE.l ) THEN
397 *
398 * RQ factorization of ( S11 S12 ) = ( 0 S12 )*Z
399 *
400  CALL cgerq2( l, n, b, ldb, tau, work, info )
401 *
402 * Update A := A*Z**H
403 *
404  CALL cunmr2( 'Right', 'Conjugate transpose', m, n, l, b, ldb,
405  $ tau, a, lda, work, info )
406  IF( wantq ) THEN
407 *
408 * Update Q := Q*Z**H
409 *
410  CALL cunmr2( 'Right', 'Conjugate transpose', n, n, l, b,
411  $ ldb, tau, q, ldq, work, info )
412  END IF
413 *
414 * Clean up B
415 *
416  CALL claset( 'Full', l, n-l, czero, czero, b, ldb )
417  DO 60 j = n - l + 1, n
418  DO 50 i = j - n + l + 1, l
419  b( i, j ) = czero
420  50 CONTINUE
421  60 CONTINUE
422 *
423  END IF
424 *
425 * Let N-L L
426 * A = ( A11 A12 ) M,
427 *
428 * then the following does the complete QR decomposition of A11:
429 *
430 * A11 = U*( 0 T12 )*P1**H
431 * ( 0 0 )
432 *
433  DO 70 i = 1, n - l
434  iwork( i ) = 0
435  70 CONTINUE
436  CALL cgeqpf( m, n-l, a, lda, iwork, tau, work, rwork, info )
437 *
438 * Determine the effective rank of A11
439 *
440  k = 0
441  DO 80 i = 1, min( m, n-l )
442  IF( cabs1( a( i, i ) ).GT.tola )
443  $ k = k + 1
444  80 CONTINUE
445 *
446 * Update A12 := U**H*A12, where A12 = A( 1:M, N-L+1:N )
447 *
448  CALL cunm2r( 'Left', 'Conjugate transpose', m, l, min( m, n-l ),
449  $ a, lda, tau, a( 1, n-l+1 ), lda, work, info )
450 *
451  IF( wantu ) THEN
452 *
453 * Copy the details of U, and form U
454 *
455  CALL claset( 'Full', m, m, czero, czero, u, ldu )
456  IF( m.GT.1 )
457  $ CALL clacpy( 'Lower', m-1, n-l, a( 2, 1 ), lda, u( 2, 1 ),
458  $ ldu )
459  CALL cung2r( m, m, min( m, n-l ), u, ldu, tau, work, info )
460  END IF
461 *
462  IF( wantq ) THEN
463 *
464 * Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1
465 *
466  CALL clapmt( forwrd, n, n-l, q, ldq, iwork )
467  END IF
468 *
469 * Clean up A: set the strictly lower triangular part of
470 * A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0.
471 *
472  DO 100 j = 1, k - 1
473  DO 90 i = j + 1, k
474  a( i, j ) = czero
475  90 CONTINUE
476  100 CONTINUE
477  IF( m.GT.k )
478  $ CALL claset( 'Full', m-k, n-l, czero, czero, a( k+1, 1 ), lda )
479 *
480  IF( n-l.GT.k ) THEN
481 *
482 * RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1
483 *
484  CALL cgerq2( k, n-l, a, lda, tau, work, info )
485 *
486  IF( wantq ) THEN
487 *
488 * Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**H
489 *
490  CALL cunmr2( 'Right', 'Conjugate transpose', n, n-l, k, a,
491  $ lda, tau, q, ldq, work, info )
492  END IF
493 *
494 * Clean up A
495 *
496  CALL claset( 'Full', k, n-l-k, czero, czero, a, lda )
497  DO 120 j = n - l - k + 1, n - l
498  DO 110 i = j - n + l + k + 1, k
499  a( i, j ) = czero
500  110 CONTINUE
501  120 CONTINUE
502 *
503  END IF
504 *
505  IF( m.GT.k ) THEN
506 *
507 * QR factorization of A( K+1:M,N-L+1:N )
508 *
509  CALL cgeqr2( m-k, l, a( k+1, n-l+1 ), lda, tau, work, info )
510 *
511  IF( wantu ) THEN
512 *
513 * Update U(:,K+1:M) := U(:,K+1:M)*U1
514 *
515  CALL cunm2r( 'Right', 'No transpose', m, m-k, min( m-k, l ),
516  $ a( k+1, n-l+1 ), lda, tau, u( 1, k+1 ), ldu,
517  $ work, info )
518  END IF
519 *
520 * Clean up
521 *
522  DO 140 j = n - l + 1, n
523  DO 130 i = j - n + k + l + 1, m
524  a( i, j ) = czero
525  130 CONTINUE
526  140 CONTINUE
527 *
528  END IF
529 *
530  RETURN
531 *
532 * End of CGGSVP
533 *
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 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 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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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 cung2r(M, N, K, A, LDA, TAU, WORK, INFO)
CUNG2R
Definition: cung2r.f:116
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
subroutine cgeqpf(M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO)
CGEQPF
Definition: cgeqpf.f:150
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

Here is the call graph for this function:

Here is the caller graph for this function: