LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zggsvp3 ( 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  LWORK,
integer  INFO 
)

ZGGSVP3

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

Purpose: 
 ZGGSVP3 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
 ZGGSVD3.
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(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 ZGEQP3 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.

ZGGSVP3 replaces the deprecated subroutine ZGGSVP.

Definition at line 280 of file zggsvp3.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  DOUBLE PRECISION tola, tolb
293 * ..
294 * .. Array Arguments ..
295  INTEGER iwork( * )
296  DOUBLE PRECISION rwork( * )
297  COMPLEX*16 a( lda, * ), b( ldb, * ), q( ldq, * ),
298  $ tau( * ), u( ldu, * ), v( ldv, * ), work( * )
299 * ..
300 *
301 * =====================================================================
302 *
303 * .. Parameters ..
304  COMPLEX*16 czero, cone
305  parameter ( czero = ( 0.0d+0, 0.0d+0 ),
306  $ cone = ( 1.0d+0, 0.0d+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 xerbla, zgeqp3, zgeqr2, zgerq2, zlacpy, zlapmt,
318  $ zlaset, zung2r, zunm2r, zunmr2
319 * ..
320 * .. Intrinsic Functions ..
321  INTRINSIC abs, dble, dimag, max, min
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 zgeqp3( 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 zgeqp3( 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 ) = dcmplx( lwkopt )
380  END IF
381 *
382  IF( info.NE.0 ) THEN
383  CALL xerbla( 'ZGGSVP3', -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 zgeqp3( p, n, b, ldb, iwork, tau, work, lwork, rwork, info )
397 *
398 * Update A := A*P
399 *
400  CALL zlapmt( 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 zlaset( 'Full', p, p, czero, czero, v, ldv )
415  IF( p.GT.1 )
416  $ CALL zlacpy( 'Lower', p-1, n, b( 2, 1 ), ldb, v( 2, 1 ),
417  $ ldv )
418  CALL zung2r( 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 zlaset( '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 zlaset( 'Full', n, n, czero, cone, q, ldq )
436  CALL zlapmt( 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 zgerq2( l, n, b, ldb, tau, work, info )
444 *
445 * Update A := A*Z**H
446 *
447  CALL zunmr2( '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 zunmr2( '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 zlaset( '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 zgeqp3( 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 zunm2r( '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 zlaset( 'Full', m, m, czero, czero, u, ldu )
500  IF( m.GT.1 )
501  $ CALL zlacpy( 'Lower', m-1, n-l, a( 2, 1 ), lda, u( 2, 1 ),
502  $ ldu )
503  CALL zung2r( 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 zlapmt( 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 zlaset( '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 zgerq2( 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 zunmr2( '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 zlaset( '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 zgeqr2( 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 zunm2r( '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 ) = dcmplx( lwkopt )
575  RETURN
576 *
577 * End of ZGGSVP3
578 *
zung2r
subroutine zung2r(M, N, K, A, LDA, TAU, WORK, INFO)
ZUNG2R
Definition: zung2r.f:116
zlacpy
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
zunmr2
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
zlaset
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
zlapmt
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
zgeqr2
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
zunm2r
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
zgeqp3
subroutine zgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
ZGEQP3
Definition: zgeqp3.f:161
zgerq2
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
lsame
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: