SRC/cggsvd3.f

   1 *> \brief <b> CGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CGGSVD3 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggsvd3.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggsvd3.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggsvd3.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CGGSVD3( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
  22 *                           LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK,
  23 *                           LWORK, RWORK, IWORK, 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 *       ..
  29 *       .. Array Arguments ..
  30 *       INTEGER            IWORK( * )
  31 *       REAL               ALPHA( * ), BETA( * ), RWORK( * )
  32 *       COMPLEX            A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
  33 *      $                   U( LDU, * ), V( LDV, * ), WORK( * )
  34 *       ..
  35 *
  36 *
  37 *> \par Purpose:
  38 *  =============
  39 *>
  40 *> \verbatim
  41 *>
  42 *> CGGSVD3 computes the generalized singular value decomposition (GSVD)
  43 *> of an M-by-N complex matrix A and P-by-N complex matrix B:
  44 *>
  45 *>       U**H*A*Q = D1*( 0 R ),    V**H*B*Q = D2*( 0 R )
  46 *>
  47 *> where U, V and Q are unitary matrices.
  48 *> Let K+L = the effective numerical rank of the
  49 *> matrix (A**H,B**H)**H, then R is a (K+L)-by-(K+L) nonsingular upper
  50 *> triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal"
  51 *> matrices and of the following structures, respectively:
  52 *>
  53 *> If M-K-L >= 0,
  54 *>
  55 *>                     K  L
  56 *>        D1 =     K ( I  0 )
  57 *>                 L ( 0  C )
  58 *>             M-K-L ( 0  0 )
  59 *>
  60 *>                   K  L
  61 *>        D2 =   L ( 0  S )
  62 *>             P-L ( 0  0 )
  63 *>
  64 *>                 N-K-L  K    L
  65 *>   ( 0 R ) = K (  0   R11  R12 )
  66 *>             L (  0    0   R22 )
  67 *>
  68 *> where
  69 *>
  70 *>   C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
  71 *>   S = diag( BETA(K+1),  ... , BETA(K+L) ),
  72 *>   C**2 + S**2 = I.
  73 *>
  74 *>   R is stored in A(1:K+L,N-K-L+1:N) on exit.
  75 *>
  76 *> If M-K-L < 0,
  77 *>
  78 *>                   K M-K K+L-M
  79 *>        D1 =   K ( I  0    0   )
  80 *>             M-K ( 0  C    0   )
  81 *>
  82 *>                     K M-K K+L-M
  83 *>        D2 =   M-K ( 0  S    0  )
  84 *>             K+L-M ( 0  0    I  )
  85 *>               P-L ( 0  0    0  )
  86 *>
  87 *>                    N-K-L  K   M-K  K+L-M
  88 *>   ( 0 R ) =     K ( 0    R11  R12  R13  )
  89 *>               M-K ( 0     0   R22  R23  )
  90 *>             K+L-M ( 0     0    0   R33  )
  91 *>
  92 *> where
  93 *>
  94 *>   C = diag( ALPHA(K+1), ... , ALPHA(M) ),
  95 *>   S = diag( BETA(K+1),  ... , BETA(M) ),
  96 *>   C**2 + S**2 = I.
  97 *>
  98 *>   (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
  99 *>   ( 0  R22 R23 )
 100 *>   in B(M-K+1:L,N+M-K-L+1:N) on exit.
 101 *>
 102 *> The routine computes C, S, R, and optionally the unitary
 103 *> transformation matrices U, V and Q.
 104 *>
 105 *> In particular, if B is an N-by-N nonsingular matrix, then the GSVD of
 106 *> A and B implicitly gives the SVD of A*inv(B):
 107 *>                      A*inv(B) = U*(D1*inv(D2))*V**H.
 108 *> If ( A**H,B**H)**H has orthonormal columns, then the GSVD of A and B is also
 109 *> equal to the CS decomposition of A and B. Furthermore, the GSVD can
 110 *> be used to derive the solution of the eigenvalue problem:
 111 *>                      A**H*A x = lambda* B**H*B x.
 112 *> In some literature, the GSVD of A and B is presented in the form
 113 *>                  U**H*A*X = ( 0 D1 ),   V**H*B*X = ( 0 D2 )
 114 *> where U and V are orthogonal and X is nonsingular, and D1 and D2 are
 115 *> ``diagonal''.  The former GSVD form can be converted to the latter
 116 *> form by taking the nonsingular matrix X as
 117 *>
 118 *>                       X = Q*(  I   0    )
 119 *>                             (  0 inv(R) )
 120 *> \endverbatim
 121 *
 122 *  Arguments:
 123 *  ==========
 124 *
 125 *> \param[in] JOBU
 126 *> \verbatim
 127 *>          JOBU is CHARACTER*1
 128 *>          = 'U':  Unitary matrix U is computed;
 129 *>          = 'N':  U is not computed.
 130 *> \endverbatim
 131 *>
 132 *> \param[in] JOBV
 133 *> \verbatim
 134 *>          JOBV is CHARACTER*1
 135 *>          = 'V':  Unitary matrix V is computed;
 136 *>          = 'N':  V is not computed.
 137 *> \endverbatim
 138 *>
 139 *> \param[in] JOBQ
 140 *> \verbatim
 141 *>          JOBQ is CHARACTER*1
 142 *>          = 'Q':  Unitary matrix Q is computed;
 143 *>          = 'N':  Q is not computed.
 144 *> \endverbatim
 145 *>
 146 *> \param[in] M
 147 *> \verbatim
 148 *>          M is INTEGER
 149 *>          The number of rows of the matrix A.  M >= 0.
 150 *> \endverbatim
 151 *>
 152 *> \param[in] N
 153 *> \verbatim
 154 *>          N is INTEGER
 155 *>          The number of columns of the matrices A and B.  N >= 0.
 156 *> \endverbatim
 157 *>
 158 *> \param[in] P
 159 *> \verbatim
 160 *>          P is INTEGER
 161 *>          The number of rows of the matrix B.  P >= 0.
 162 *> \endverbatim
 163 *>
 164 *> \param[out] K
 165 *> \verbatim
 166 *>          K is INTEGER
 167 *> \endverbatim
 168 *>
 169 *> \param[out] L
 170 *> \verbatim
 171 *>          L is INTEGER
 172 *>
 173 *>          On exit, K and L specify the dimension of the subblocks
 174 *>          described in Purpose.
 175 *>          K + L = effective numerical rank of (A**H,B**H)**H.
 176 *> \endverbatim
 177 *>
 178 *> \param[in,out] A
 179 *> \verbatim
 180 *>          A is COMPLEX array, dimension (LDA,N)
 181 *>          On entry, the M-by-N matrix A.
 182 *>          On exit, A contains the triangular matrix R, or part of R.
 183 *>          See Purpose for details.
 184 *> \endverbatim
 185 *>
 186 *> \param[in] LDA
 187 *> \verbatim
 188 *>          LDA is INTEGER
 189 *>          The leading dimension of the array A. LDA >= max(1,M).
 190 *> \endverbatim
 191 *>
 192 *> \param[in,out] B
 193 *> \verbatim
 194 *>          B is COMPLEX array, dimension (LDB,N)
 195 *>          On entry, the P-by-N matrix B.
 196 *>          On exit, B contains part of the triangular matrix R if
 197 *>          M-K-L < 0.  See Purpose for details.
 198 *> \endverbatim
 199 *>
 200 *> \param[in] LDB
 201 *> \verbatim
 202 *>          LDB is INTEGER
 203 *>          The leading dimension of the array B. LDB >= max(1,P).
 204 *> \endverbatim
 205 *>
 206 *> \param[out] ALPHA
 207 *> \verbatim
 208 *>          ALPHA is REAL array, dimension (N)
 209 *> \endverbatim
 210 *>
 211 *> \param[out] BETA
 212 *> \verbatim
 213 *>          BETA is REAL array, dimension (N)
 214 *>
 215 *>          On exit, ALPHA and BETA contain the generalized singular
 216 *>          value pairs of A and B;
 217 *>            ALPHA(1:K) = 1,
 218 *>            BETA(1:K)  = 0,
 219 *>          and if M-K-L >= 0,
 220 *>            ALPHA(K+1:K+L) = C,
 221 *>            BETA(K+1:K+L)  = S,
 222 *>          or if M-K-L < 0,
 223 *>            ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0
 224 *>            BETA(K+1:M) =S, BETA(M+1:K+L) =1
 225 *>          and
 226 *>            ALPHA(K+L+1:N) = 0
 227 *>            BETA(K+L+1:N)  = 0
 228 *> \endverbatim
 229 *>
 230 *> \param[out] U
 231 *> \verbatim
 232 *>          U is COMPLEX array, dimension (LDU,M)
 233 *>          If JOBU = 'U', U contains the M-by-M unitary matrix U.
 234 *>          If JOBU = 'N', U is not referenced.
 235 *> \endverbatim
 236 *>
 237 *> \param[in] LDU
 238 *> \verbatim
 239 *>          LDU is INTEGER
 240 *>          The leading dimension of the array U. LDU >= max(1,M) if
 241 *>          JOBU = 'U'; LDU >= 1 otherwise.
 242 *> \endverbatim
 243 *>
 244 *> \param[out] V
 245 *> \verbatim
 246 *>          V is COMPLEX array, dimension (LDV,P)
 247 *>          If JOBV = 'V', V contains the P-by-P unitary matrix V.
 248 *>          If JOBV = 'N', V is not referenced.
 249 *> \endverbatim
 250 *>
 251 *> \param[in] LDV
 252 *> \verbatim
 253 *>          LDV is INTEGER
 254 *>          The leading dimension of the array V. LDV >= max(1,P) if
 255 *>          JOBV = 'V'; LDV >= 1 otherwise.
 256 *> \endverbatim
 257 *>
 258 *> \param[out] Q
 259 *> \verbatim
 260 *>          Q is COMPLEX array, dimension (LDQ,N)
 261 *>          If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q.
 262 *>          If JOBQ = 'N', Q is not referenced.
 263 *> \endverbatim
 264 *>
 265 *> \param[in] LDQ
 266 *> \verbatim
 267 *>          LDQ is INTEGER
 268 *>          The leading dimension of the array Q. LDQ >= max(1,N) if
 269 *>          JOBQ = 'Q'; LDQ >= 1 otherwise.
 270 *> \endverbatim
 271 *>
 272 *> \param[out] WORK
 273 *> \verbatim
 274 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
 275 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 276 *> \endverbatim
 277 *>
 278 *> \param[in] LWORK
 279 *> \verbatim
 280 *>          LWORK is INTEGER
 281 *>          The dimension of the array WORK.
 282 *>
 283 *>          If LWORK = -1, then a workspace query is assumed; the routine
 284 *>          only calculates the optimal size of the WORK array, returns
 285 *>          this value as the first entry of the WORK array, and no error
 286 *>          message related to LWORK is issued by XERBLA.
 287 *> \endverbatim
 288 *>
 289 *> \param[out] RWORK
 290 *> \verbatim
 291 *>          RWORK is REAL array, dimension (2*N)
 292 *> \endverbatim
 293 *>
 294 *> \param[out] IWORK
 295 *> \verbatim
 296 *>          IWORK is INTEGER array, dimension (N)
 297 *>          On exit, IWORK stores the sorting information. More
 298 *>          precisely, the following loop will sort ALPHA
 299 *>             for I = K+1, min(M,K+L)
 300 *>                 swap ALPHA(I) and ALPHA(IWORK(I))
 301 *>             endfor
 302 *>          such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N).
 303 *> \endverbatim
 304 *>
 305 *> \param[out] INFO
 306 *> \verbatim
 307 *>          INFO is INTEGER
 308 *>          = 0:  successful exit.
 309 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 310 *>          > 0:  if INFO = 1, the Jacobi-type procedure failed to
 311 *>                converge.  For further details, see subroutine CTGSJA.
 312 *> \endverbatim
 313 *
 314 *> \par Internal Parameters:
 315 *  =========================
 316 *>
 317 *> \verbatim
 318 *>  TOLA    REAL
 319 *>  TOLB    REAL
 320 *>          TOLA and TOLB are the thresholds to determine the effective
 321 *>          rank of (A**H,B**H)**H. Generally, they are set to
 322 *>                   TOLA = MAX(M,N)*norm(A)*MACHEPS,
 323 *>                   TOLB = MAX(P,N)*norm(B)*MACHEPS.
 324 *>          The size of TOLA and TOLB may affect the size of backward
 325 *>          errors of the decomposition.
 326 *> \endverbatim
 327 *
 328 *  Authors:
 329 *  ========
 330 *
 331 *> \author Univ. of Tennessee
 332 *> \author Univ. of California Berkeley
 333 *> \author Univ. of Colorado Denver
 334 *> \author NAG Ltd.
 335 *
 336 *> \date August 2015
 337 *
 338 *> \ingroup complexOTHERsing
 339 *
 340 *> \par Contributors:
 341 *  ==================
 342 *>
 343 *>     Ming Gu and Huan Ren, Computer Science Division, University of
 344 *>     California at Berkeley, USA
 345 *>
 346 *
 347 *> \par Further Details:
 348 *  =====================
 349 *>
 350 *>  CGGSVD3 replaces the deprecated subroutine CGGSVD.
 351 *>
 352 *  =====================================================================
 353       SUBROUTINE CGGSVD3( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
 354      $                    LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ,
 355      $                    WORK, LWORK, RWORK, IWORK, INFO )
 356 *
 357 *  -- LAPACK driver routine (version 3.6.0) --
 358 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 359 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 360 *     August 2015
 361 *
 362 *     .. Scalar Arguments ..
 363       CHARACTER          JOBQ, JOBU, JOBV
 364       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P,
 365      $                   LWORK
 366 *     ..
 367 *     .. Array Arguments ..
 368       INTEGER            IWORK( * )
 369       REAL               ALPHA( * ), BETA( * ), RWORK( * )
 370       COMPLEX            A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
 371      $                   U( LDU, * ), V( LDV, * ), WORK( * )
 372 *     ..
 373 *
 374 *  =====================================================================
 375 *
 376 *     .. Local Scalars ..
 377       LOGICAL            WANTQ, WANTU, WANTV, LQUERY
 378       INTEGER            I, IBND, ISUB, J, NCYCLE, LWKOPT
 379       REAL               ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
 380 *     ..
 381 *     .. External Functions ..
 382       LOGICAL            LSAME
 383       REAL               CLANGE, SLAMCH
 384       EXTERNAL           LSAME, CLANGE, SLAMCH
 385 *     ..
 386 *     .. External Subroutines ..
 387       EXTERNAL           CGGSVP3, CTGSJA, SCOPY, XERBLA
 388 *     ..
 389 *     .. Intrinsic Functions ..
 390       INTRINSIC          MAX, MIN
 391 *     ..
 392 *     .. Executable Statements ..
 393 *
 394 *     Decode and test the input parameters
 395 *
 396       WANTU = LSAME( JOBU, 'U' )
 397       WANTV = LSAME( JOBV, 'V' )
 398       WANTQ = LSAME( JOBQ, 'Q' )
 399       LQUERY = ( LWORK.EQ.-1 )
 400       LWKOPT = 1
 401 *
 402 *     Test the input arguments
 403 *
 404       INFO = 0
 405       IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN
 406          INFO = -1
 407       ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN
 408          INFO = -2
 409       ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN
 410          INFO = -3
 411       ELSE IF( M.LT.0 ) THEN
 412          INFO = -4
 413       ELSE IF( N.LT.0 ) THEN
 414          INFO = -5
 415       ELSE IF( P.LT.0 ) THEN
 416          INFO = -6
 417       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 418          INFO = -10
 419       ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
 420          INFO = -12
 421       ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
 422          INFO = -16
 423       ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
 424          INFO = -18
 425       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
 426          INFO = -20
 427       ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
 428          INFO = -24
 429       END IF
 430 *
 431 *     Compute workspace
 432 *
 433       IF( INFO.EQ.0 ) THEN
 434          CALL CGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA,
 435      $                 TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK,
 436      $                 WORK, WORK, -1, INFO )
 437          LWKOPT = N + INT( WORK( 1 ) )
 438          LWKOPT = MAX( 2*N, LWKOPT )
 439          LWKOPT = MAX( 1, LWKOPT )
 440          WORK( 1 ) = CMPLX( LWKOPT )
 441       END IF
 442 *
 443       IF( INFO.NE.0 ) THEN
 444          CALL XERBLA( 'CGGSVD3', -INFO )
 445          RETURN
 446       END IF
 447       IF( LQUERY ) THEN
 448          RETURN
 449       ENDIF
 450 *
 451 *     Compute the Frobenius norm of matrices A and B
 452 *
 453       ANORM = CLANGE( '1', M, N, A, LDA, RWORK )
 454       BNORM = CLANGE( '1', P, N, B, LDB, RWORK )
 455 *
 456 *     Get machine precision and set up threshold for determining
 457 *     the effective numerical rank of the matrices A and B.
 458 *
 459       ULP = SLAMCH( 'Precision' )
 460       UNFL = SLAMCH( 'Safe Minimum' )
 461       TOLA = MAX( M, N )*MAX( ANORM, UNFL )*ULP
 462       TOLB = MAX( P, N )*MAX( BNORM, UNFL )*ULP
 463 *
 464       CALL CGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA,
 465      $              TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK,
 466      $              WORK, WORK( N+1 ), LWORK-N, INFO )
 467 *
 468 *     Compute the GSVD of two upper "triangular" matrices
 469 *
 470       CALL CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB,
 471      $             TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ,
 472      $             WORK, NCYCLE, INFO )
 473 *
 474 *     Sort the singular values and store the pivot indices in IWORK
 475 *     Copy ALPHA to RWORK, then sort ALPHA in RWORK
 476 *
 477       CALL SCOPY( N, ALPHA, 1, RWORK, 1 )
 478       IBND = MIN( L, M-K )
 479       DO 20 I = 1, IBND
 480 *
 481 *        Scan for largest ALPHA(K+I)
 482 *
 483          ISUB = I
 484          SMAX = RWORK( K+I )
 485          DO 10 J = I + 1, IBND
 486             TEMP = RWORK( K+J )
 487             IF( TEMP.GT.SMAX ) THEN
 488                ISUB = J
 489                SMAX = TEMP
 490             END IF
 491    10    CONTINUE
 492          IF( ISUB.NE.I ) THEN
 493             RWORK( K+ISUB ) = RWORK( K+I )
 494             RWORK( K+I ) = SMAX
 495             IWORK( K+I ) = K + ISUB
 496          ELSE
 497             IWORK( K+I ) = K + I
 498          END IF
 499    20 CONTINUE
 500 *
 501       WORK( 1 ) = CMPLX( LWKOPT )
 502       RETURN
 503 *
 504 *     End of CGGSVD3
 505 *
 506       END