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

Go to the source code of this file.

Functions/Subroutines

subroutine zggsvp (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)
 ZGGSVP More...
 

Function/Subroutine Documentation

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

ZGGSVP

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

Purpose:
 ZGGSVP 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
 ZGGSVD.
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*16 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*16 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 DOUBLE PRECISION
[in]TOLB
          TOLB is DOUBLE PRECISION

          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)*MAZHEPS,
             TOLB = MAX(P,N)*norm(B)*MAZHEPS.
          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*16 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*16 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*16 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 DOUBLE PRECISION array, dimension (2*N)
[out]TAU
          TAU is COMPLEX*16 array, dimension (N)
[out]WORK
          WORK is COMPLEX*16 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 ZGEQPF 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 265 of file zggsvp.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: