LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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,
315 $ zlaset, zung2r, zunm2r, zunmr2
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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zgeqp3
subroutine zgeqp3(m, n, a, lda, jpvt, tau, work, lwork, rwork, info)
ZGEQP3
Definition zgeqp3.f:159
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:130
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:123
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:103
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:104
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:106
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
zung2r
subroutine zung2r(m, n, k, a, lda, tau, work, info)
ZUNG2R
Definition zung2r.f:114
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:159
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:159
Here is the call graph for this function:
Here is the caller graph for this function: