LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zggsvp3()

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.
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 275 of file zggsvp3.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  DOUBLE PRECISION TOLA, TOLB
290 * ..
291 * .. Array Arguments ..
292  INTEGER IWORK( * )
293  DOUBLE PRECISION RWORK( * )
294  COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
295  $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
296 * ..
297 *
298 * =====================================================================
299 *
300 * .. Parameters ..
301  COMPLEX*16 CZERO, CONE
302  parameter( czero = ( 0.0d+0, 0.0d+0 ),
303  $ cone = ( 1.0d+0, 0.0d+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 xerbla, zgeqp3, zgeqr2, zgerq2, zlacpy, zlapmt,
316 * ..
317 * .. Intrinsic Functions ..
318  INTRINSIC abs, dble, dimag, max, min
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 zgeqp3( 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 zgeqp3( 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 ) = dcmplx( lwkopt )
377  END IF
378 *
379  IF( info.NE.0 ) THEN
380  CALL xerbla( 'ZGGSVP3', -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 zgeqp3( p, n, b, ldb, iwork, tau, work, lwork, rwork, info )
394 *
395 * Update A := A*P
396 *
397  CALL zlapmt( 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 zlaset( 'Full', p, p, czero, czero, v, ldv )
412  IF( p.GT.1 )
413  $ CALL zlacpy( 'Lower', p-1, n, b( 2, 1 ), ldb, v( 2, 1 ),
414  $ ldv )
415  CALL zung2r( 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 zlaset( '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 zlaset( 'Full', n, n, czero, cone, q, ldq )
433  CALL zlapmt( 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 zgerq2( l, n, b, ldb, tau, work, info )
441 *
442 * Update A := A*Z**H
443 *
444  CALL zunmr2( '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 zunmr2( '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 zlaset( '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 zgeqp3( 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 zunm2r( '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 zlaset( 'Full', m, m, czero, czero, u, ldu )
497  IF( m.GT.1 )
498  $ CALL zlacpy( 'Lower', m-1, n-l, a( 2, 1 ), lda, u( 2, 1 ),
499  $ ldu )
500  CALL zung2r( 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 zlapmt( 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 zlaset( '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 zgerq2( 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 zunmr2( '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 zlaset( '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 zgeqr2( 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 zunm2r( '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 ) = dcmplx( lwkopt )
572  RETURN
573 *
574 * End of ZGGSVP3
575 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
ZGEQP3
Definition: zgeqp3.f:159
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:130
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:123
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zlapmt(FORWRD, M, N, X, LDX, K)
ZLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: zlapmt.f:104
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:106
subroutine zung2r(M, N, K, A, LDA, TAU, WORK, INFO)
ZUNG2R
Definition: zung2r.f:114
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:159
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:159
Here is the call graph for this function:
Here is the caller graph for this function: