LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cggsvp3()

subroutine cggsvp3 ( 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  LWORK,
integer  INFO 
)

CGGSVP3

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

Purpose:
 CGGSVP3 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
 CGGSVD3.
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(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 CGEQP3 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.

  CGGSVP3 replaces the deprecated subroutine CGGSVP.

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