SRC/cggsvp3.f

   1 *> \brief \b CGGSVP3
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CGGSVP3 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggsvp3.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggsvp3.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggsvp3.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
  22 *                           TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
  23 *                           IWORK, RWORK, TAU, WORK, LWORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          JOBQ, JOBU, JOBV
  27 *       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, LWORK
  28 *       REAL               TOLA, TOLB
  29 *       ..
  30 *       .. Array Arguments ..
  31 *       INTEGER            IWORK( * )
  32 *       REAL               RWORK( * )
  33 *       COMPLEX            A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
  34 *      $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
  35 *       ..
  36 *
  37 *
  38 *> \par Purpose:
  39 *  =============
  40 *>
  41 *> \verbatim
  42 *>
  43 *> CGGSVP3 computes unitary matrices U, V and Q such that
  44 *>
  45 *>                    N-K-L  K    L
  46 *>  U**H*A*Q =     K ( 0    A12  A13 )  if M-K-L >= 0;
  47 *>                 L ( 0     0   A23 )
  48 *>             M-K-L ( 0     0    0  )
  49 *>
  50 *>                  N-K-L  K    L
  51 *>         =     K ( 0    A12  A13 )  if M-K-L < 0;
  52 *>             M-K ( 0     0   A23 )
  53 *>
  54 *>                  N-K-L  K    L
  55 *>  V**H*B*Q =   L ( 0     0   B13 )
  56 *>             P-L ( 0     0    0  )
  57 *>
  58 *> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
  59 *> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
  60 *> otherwise A23 is (M-K)-by-L upper trapezoidal.  K+L = the effective
  61 *> numerical rank of the (M+P)-by-N matrix (A**H,B**H)**H.
  62 *>
  63 *> This decomposition is the preprocessing step for computing the
  64 *> Generalized Singular Value Decomposition (GSVD), see subroutine
  65 *> CGGSVD3.
  66 *> \endverbatim
  67 *
  68 *  Arguments:
  69 *  ==========
  70 *
  71 *> \param[in] JOBU
  72 *> \verbatim
  73 *>          JOBU is CHARACTER*1
  74 *>          = 'U':  Unitary matrix U is computed;
  75 *>          = 'N':  U is not computed.
  76 *> \endverbatim
  77 *>
  78 *> \param[in] JOBV
  79 *> \verbatim
  80 *>          JOBV is CHARACTER*1
  81 *>          = 'V':  Unitary matrix V is computed;
  82 *>          = 'N':  V is not computed.
  83 *> \endverbatim
  84 *>
  85 *> \param[in] JOBQ
  86 *> \verbatim
  87 *>          JOBQ is CHARACTER*1
  88 *>          = 'Q':  Unitary matrix Q is computed;
  89 *>          = 'N':  Q is not computed.
  90 *> \endverbatim
  91 *>
  92 *> \param[in] M
  93 *> \verbatim
  94 *>          M is INTEGER
  95 *>          The number of rows of the matrix A.  M >= 0.
  96 *> \endverbatim
  97 *>
  98 *> \param[in] P
  99 *> \verbatim
 100 *>          P is INTEGER
 101 *>          The number of rows of the matrix B.  P >= 0.
 102 *> \endverbatim
 103 *>
 104 *> \param[in] N
 105 *> \verbatim
 106 *>          N is INTEGER
 107 *>          The number of columns of the matrices A and B.  N >= 0.
 108 *> \endverbatim
 109 *>
 110 *> \param[in,out] A
 111 *> \verbatim
 112 *>          A is COMPLEX array, dimension (LDA,N)
 113 *>          On entry, the M-by-N matrix A.
 114 *>          On exit, A contains the triangular (or trapezoidal) matrix
 115 *>          described in the Purpose section.
 116 *> \endverbatim
 117 *>
 118 *> \param[in] LDA
 119 *> \verbatim
 120 *>          LDA is INTEGER
 121 *>          The leading dimension of the array A. LDA >= max(1,M).
 122 *> \endverbatim
 123 *>
 124 *> \param[in,out] B
 125 *> \verbatim
 126 *>          B is COMPLEX array, dimension (LDB,N)
 127 *>          On entry, the P-by-N matrix B.
 128 *>          On exit, B contains the triangular matrix described in
 129 *>          the Purpose section.
 130 *> \endverbatim
 131 *>
 132 *> \param[in] LDB
 133 *> \verbatim
 134 *>          LDB is INTEGER
 135 *>          The leading dimension of the array B. LDB >= max(1,P).
 136 *> \endverbatim
 137 *>
 138 *> \param[in] TOLA
 139 *> \verbatim
 140 *>          TOLA is REAL
 141 *> \endverbatim
 142 *>
 143 *> \param[in] TOLB
 144 *> \verbatim
 145 *>          TOLB is REAL
 146 *>
 147 *>          TOLA and TOLB are the thresholds to determine the effective
 148 *>          numerical rank of matrix B and a subblock of A. Generally,
 149 *>          they are set to
 150 *>             TOLA = MAX(M,N)*norm(A)*MACHEPS,
 151 *>             TOLB = MAX(P,N)*norm(B)*MACHEPS.
 152 *>          The size of TOLA and TOLB may affect the size of backward
 153 *>          errors of the decomposition.
 154 *> \endverbatim
 155 *>
 156 *> \param[out] K
 157 *> \verbatim
 158 *>          K is INTEGER
 159 *> \endverbatim
 160 *>
 161 *> \param[out] L
 162 *> \verbatim
 163 *>          L is INTEGER
 164 *>
 165 *>          On exit, K and L specify the dimension of the subblocks
 166 *>          described in Purpose section.
 167 *>          K + L = effective numerical rank of (A**H,B**H)**H.
 168 *> \endverbatim
 169 *>
 170 *> \param[out] U
 171 *> \verbatim
 172 *>          U is COMPLEX array, dimension (LDU,M)
 173 *>          If JOBU = 'U', U contains the unitary matrix U.
 174 *>          If JOBU = 'N', U is not referenced.
 175 *> \endverbatim
 176 *>
 177 *> \param[in] LDU
 178 *> \verbatim
 179 *>          LDU is INTEGER
 180 *>          The leading dimension of the array U. LDU >= max(1,M) if
 181 *>          JOBU = 'U'; LDU >= 1 otherwise.
 182 *> \endverbatim
 183 *>
 184 *> \param[out] V
 185 *> \verbatim
 186 *>          V is COMPLEX array, dimension (LDV,P)
 187 *>          If JOBV = 'V', V contains the unitary matrix V.
 188 *>          If JOBV = 'N', V is not referenced.
 189 *> \endverbatim
 190 *>
 191 *> \param[in] LDV
 192 *> \verbatim
 193 *>          LDV is INTEGER
 194 *>          The leading dimension of the array V. LDV >= max(1,P) if
 195 *>          JOBV = 'V'; LDV >= 1 otherwise.
 196 *> \endverbatim
 197 *>
 198 *> \param[out] Q
 199 *> \verbatim
 200 *>          Q is COMPLEX array, dimension (LDQ,N)
 201 *>          If JOBQ = 'Q', Q contains the unitary matrix Q.
 202 *>          If JOBQ = 'N', Q is not referenced.
 203 *> \endverbatim
 204 *>
 205 *> \param[in] LDQ
 206 *> \verbatim
 207 *>          LDQ is INTEGER
 208 *>          The leading dimension of the array Q. LDQ >= max(1,N) if
 209 *>          JOBQ = 'Q'; LDQ >= 1 otherwise.
 210 *> \endverbatim
 211 *>
 212 *> \param[out] IWORK
 213 *> \verbatim
 214 *>          IWORK is INTEGER array, dimension (N)
 215 *> \endverbatim
 216 *>
 217 *> \param[out] RWORK
 218 *> \verbatim
 219 *>          RWORK is REAL array, dimension (2*N)
 220 *> \endverbatim
 221 *>
 222 *> \param[out] TAU
 223 *> \verbatim
 224 *>          TAU is COMPLEX array, dimension (N)
 225 *> \endverbatim
 226 *>
 227 *> \param[out] WORK
 228 *> \verbatim
 229 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
 230 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 231 *> \endverbatim
 232 *>
 233 *> \param[in] LWORK
 234 *> \verbatim
 235 *>          LWORK is INTEGER
 236 *>          The dimension of the array WORK.
 237 *>
 238 *>          If LWORK = -1, then a workspace query is assumed; the routine
 239 *>          only calculates the optimal size of the WORK array, returns
 240 *>          this value as the first entry of the WORK array, and no error
 241 *>          message related to LWORK is issued by XERBLA.
 242 *> \endverbatim
 243 *>
 244 *> \param[out] INFO
 245 *> \verbatim
 246 *>          INFO is INTEGER
 247 *>          = 0:  successful exit
 248 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 249 *> \endverbatim
 250 *
 251 *  Authors:
 252 *  ========
 253 *
 254 *> \author Univ. of Tennessee
 255 *> \author Univ. of California Berkeley
 256 *> \author Univ. of Colorado Denver
 257 *> \author NAG Ltd.
 258 *
 259 *> \date August 2015
 260 *
 261 *> \ingroup complexOTHERcomputational
 262 *
 263 *> \par Further Details:
 264 *  =====================
 265 *>
 266 *> \verbatim
 267 *>
 268 *>  The subroutine uses LAPACK subroutine CGEQP3 for the QR factorization
 269 *>  with column pivoting to detect the effective numerical rank of the
 270 *>  a matrix. It may be replaced by a better rank determination strategy.
 271 *>
 272 *>  CGGSVP3 replaces the deprecated subroutine CGGSVP.
 273 *>
 274 *> \endverbatim
 275 *>
 276 *  =====================================================================
 277       SUBROUTINE CGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
 278      $                    TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
 279      $                    IWORK, RWORK, TAU, WORK, LWORK, INFO )
 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       REAL               TOLA, TOLB
 293 *     ..
 294 *     .. Array Arguments ..
 295       INTEGER            IWORK( * )
 296       REAL               RWORK( * )
 297       COMPLEX            A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
 298      $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
 299 *     ..
 300 *
 301 *  =====================================================================
 302 *
 303 *     .. Parameters ..
 304       COMPLEX            CZERO, CONE
 305       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
 306      $                   CONE = ( 1.0E+0, 0.0E+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           CGEQP3, CGEQR2, CGERQ2, CLACPY, CLAPMT,
 318      $                   CLASET, CUNG2R, CUNM2R, CUNMR2, XERBLA
 319 *     ..
 320 *     .. Intrinsic Functions ..
 321       INTRINSIC          ABS, AIMAG, MAX, MIN, REAL
 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 CGEQP3( 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 CGEQP3( 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 ) = CMPLX( LWKOPT )
 380       END IF
 381 *
 382       IF( INFO.NE.0 ) THEN
 383          CALL XERBLA( 'CGGSVP3', -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 CGEQP3( P, N, B, LDB, IWORK, TAU, WORK, LWORK, RWORK, INFO )
 397 *
 398 *     Update A := A*P
 399 *
 400       CALL CLAPMT( 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 CLASET( 'Full', P, P, CZERO, CZERO, V, LDV )
 415          IF( P.GT.1 )
 416      $      CALL CLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ),
 417      $                   LDV )
 418          CALL CUNG2R( 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 CLASET( '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 CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
 436          CALL CLAPMT( 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 CGERQ2( L, N, B, LDB, TAU, WORK, INFO )
 444 *
 445 *        Update A := A*Z**H
 446 *
 447          CALL CUNMR2( '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 CUNMR2( '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 CLASET( '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 CGEQP3( 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 CUNM2R( '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 CLASET( 'Full', M, M, CZERO, CZERO, U, LDU )
 500          IF( M.GT.1 )
 501      $      CALL CLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ),
 502      $                   LDU )
 503          CALL CUNG2R( 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 CLAPMT( 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 CLASET( '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 CGERQ2( 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 CUNMR2( '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 CLASET( '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 CGEQR2( 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 CUNM2R( '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 ) = CMPLX( LWKOPT )
 575       RETURN
 576 *
 577 *     End of CGGSVP3
 578 *
 579       END