LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cggsvd.f
Go to the documentation of this file.
00001 *> \brief <b> CGGSVD computes the singular value decomposition (SVD) for OTHER matrices</b>
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CGGSVD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggsvd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggsvd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggsvd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
00022 *                          LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK,
00023 *                          RWORK, IWORK, INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       CHARACTER          JOBQ, JOBU, JOBV
00027 *       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       INTEGER            IWORK( * )
00031 *       REAL               ALPHA( * ), BETA( * ), RWORK( * )
00032 *       COMPLEX            A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
00033 *      $                   U( LDU, * ), V( LDV, * ), WORK( * )
00034 *       ..
00035 *  
00036 *
00037 *> \par Purpose:
00038 *  =============
00039 *>
00040 *> \verbatim
00041 *>
00042 *> CGGSVD computes the generalized singular value decomposition (GSVD)
00043 *> of an M-by-N complex matrix A and P-by-N complex matrix B:
00044 *>
00045 *>       U**H*A*Q = D1*( 0 R ),    V**H*B*Q = D2*( 0 R )
00046 *>
00047 *> where U, V and Q are unitary matrices.
00048 *> Let K+L = the effective numerical rank of the
00049 *> matrix (A**H,B**H)**H, then R is a (K+L)-by-(K+L) nonsingular upper
00050 *> triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal"
00051 *> matrices and of the following structures, respectively:
00052 *>
00053 *> If M-K-L >= 0,
00054 *>
00055 *>                     K  L
00056 *>        D1 =     K ( I  0 )
00057 *>                 L ( 0  C )
00058 *>             M-K-L ( 0  0 )
00059 *>
00060 *>                   K  L
00061 *>        D2 =   L ( 0  S )
00062 *>             P-L ( 0  0 )
00063 *>
00064 *>                 N-K-L  K    L
00065 *>   ( 0 R ) = K (  0   R11  R12 )
00066 *>             L (  0    0   R22 )
00067 *>
00068 *> where
00069 *>
00070 *>   C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
00071 *>   S = diag( BETA(K+1),  ... , BETA(K+L) ),
00072 *>   C**2 + S**2 = I.
00073 *>
00074 *>   R is stored in A(1:K+L,N-K-L+1:N) on exit.
00075 *>
00076 *> If M-K-L < 0,
00077 *>
00078 *>                   K M-K K+L-M
00079 *>        D1 =   K ( I  0    0   )
00080 *>             M-K ( 0  C    0   )
00081 *>
00082 *>                     K M-K K+L-M
00083 *>        D2 =   M-K ( 0  S    0  )
00084 *>             K+L-M ( 0  0    I  )
00085 *>               P-L ( 0  0    0  )
00086 *>
00087 *>                    N-K-L  K   M-K  K+L-M
00088 *>   ( 0 R ) =     K ( 0    R11  R12  R13  )
00089 *>               M-K ( 0     0   R22  R23  )
00090 *>             K+L-M ( 0     0    0   R33  )
00091 *>
00092 *> where
00093 *>
00094 *>   C = diag( ALPHA(K+1), ... , ALPHA(M) ),
00095 *>   S = diag( BETA(K+1),  ... , BETA(M) ),
00096 *>   C**2 + S**2 = I.
00097 *>
00098 *>   (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
00099 *>   ( 0  R22 R23 )
00100 *>   in B(M-K+1:L,N+M-K-L+1:N) on exit.
00101 *>
00102 *> The routine computes C, S, R, and optionally the unitary
00103 *> transformation matrices U, V and Q.
00104 *>
00105 *> In particular, if B is an N-by-N nonsingular matrix, then the GSVD of
00106 *> A and B implicitly gives the SVD of A*inv(B):
00107 *>                      A*inv(B) = U*(D1*inv(D2))*V**H.
00108 *> If ( A**H,B**H)**H has orthnormal columns, then the GSVD of A and B is also
00109 *> equal to the CS decomposition of A and B. Furthermore, the GSVD can
00110 *> be used to derive the solution of the eigenvalue problem:
00111 *>                      A**H*A x = lambda* B**H*B x.
00112 *> In some literature, the GSVD of A and B is presented in the form
00113 *>                  U**H*A*X = ( 0 D1 ),   V**H*B*X = ( 0 D2 )
00114 *> where U and V are orthogonal and X is nonsingular, and D1 and D2 are
00115 *> ``diagonal''.  The former GSVD form can be converted to the latter
00116 *> form by taking the nonsingular matrix X as
00117 *>
00118 *>                       X = Q*(  I   0    )
00119 *>                             (  0 inv(R) )
00120 *> \endverbatim
00121 *
00122 *  Arguments:
00123 *  ==========
00124 *
00125 *> \param[in] JOBU
00126 *> \verbatim
00127 *>          JOBU is CHARACTER*1
00128 *>          = 'U':  Unitary matrix U is computed;
00129 *>          = 'N':  U is not computed.
00130 *> \endverbatim
00131 *>
00132 *> \param[in] JOBV
00133 *> \verbatim
00134 *>          JOBV is CHARACTER*1
00135 *>          = 'V':  Unitary matrix V is computed;
00136 *>          = 'N':  V is not computed.
00137 *> \endverbatim
00138 *>
00139 *> \param[in] JOBQ
00140 *> \verbatim
00141 *>          JOBQ is CHARACTER*1
00142 *>          = 'Q':  Unitary matrix Q is computed;
00143 *>          = 'N':  Q is not computed.
00144 *> \endverbatim
00145 *>
00146 *> \param[in] M
00147 *> \verbatim
00148 *>          M is INTEGER
00149 *>          The number of rows of the matrix A.  M >= 0.
00150 *> \endverbatim
00151 *>
00152 *> \param[in] N
00153 *> \verbatim
00154 *>          N is INTEGER
00155 *>          The number of columns of the matrices A and B.  N >= 0.
00156 *> \endverbatim
00157 *>
00158 *> \param[in] P
00159 *> \verbatim
00160 *>          P is INTEGER
00161 *>          The number of rows of the matrix B.  P >= 0.
00162 *> \endverbatim
00163 *>
00164 *> \param[out] K
00165 *> \verbatim
00166 *>          K is INTEGER
00167 *> \endverbatim
00168 *>
00169 *> \param[out] L
00170 *> \verbatim
00171 *>          L is INTEGER
00172 *>
00173 *>          On exit, K and L specify the dimension of the subblocks
00174 *>          described in Purpose.
00175 *>          K + L = effective numerical rank of (A**H,B**H)**H.
00176 *> \endverbatim
00177 *>
00178 *> \param[in,out] A
00179 *> \verbatim
00180 *>          A is COMPLEX array, dimension (LDA,N)
00181 *>          On entry, the M-by-N matrix A.
00182 *>          On exit, A contains the triangular matrix R, or part of R.
00183 *>          See Purpose for details.
00184 *> \endverbatim
00185 *>
00186 *> \param[in] LDA
00187 *> \verbatim
00188 *>          LDA is INTEGER
00189 *>          The leading dimension of the array A. LDA >= max(1,M).
00190 *> \endverbatim
00191 *>
00192 *> \param[in,out] B
00193 *> \verbatim
00194 *>          B is COMPLEX array, dimension (LDB,N)
00195 *>          On entry, the P-by-N matrix B.
00196 *>          On exit, B contains part of the triangular matrix R if
00197 *>          M-K-L < 0.  See Purpose for details.
00198 *> \endverbatim
00199 *>
00200 *> \param[in] LDB
00201 *> \verbatim
00202 *>          LDB is INTEGER
00203 *>          The leading dimension of the array B. LDB >= max(1,P).
00204 *> \endverbatim
00205 *>
00206 *> \param[out] ALPHA
00207 *> \verbatim
00208 *>          ALPHA is REAL array, dimension (N)
00209 *> \endverbatim
00210 *>
00211 *> \param[out] BETA
00212 *> \verbatim
00213 *>          BETA is REAL array, dimension (N)
00214 *>
00215 *>          On exit, ALPHA and BETA contain the generalized singular
00216 *>          value pairs of A and B;
00217 *>            ALPHA(1:K) = 1,
00218 *>            BETA(1:K)  = 0,
00219 *>          and if M-K-L >= 0,
00220 *>            ALPHA(K+1:K+L) = C,
00221 *>            BETA(K+1:K+L)  = S,
00222 *>          or if M-K-L < 0,
00223 *>            ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0
00224 *>            BETA(K+1:M) =S, BETA(M+1:K+L) =1
00225 *>          and
00226 *>            ALPHA(K+L+1:N) = 0
00227 *>            BETA(K+L+1:N)  = 0
00228 *> \endverbatim
00229 *>
00230 *> \param[out] U
00231 *> \verbatim
00232 *>          U is COMPLEX array, dimension (LDU,M)
00233 *>          If JOBU = 'U', U contains the M-by-M unitary matrix U.
00234 *>          If JOBU = 'N', U is not referenced.
00235 *> \endverbatim
00236 *>
00237 *> \param[in] LDU
00238 *> \verbatim
00239 *>          LDU is INTEGER
00240 *>          The leading dimension of the array U. LDU >= max(1,M) if
00241 *>          JOBU = 'U'; LDU >= 1 otherwise.
00242 *> \endverbatim
00243 *>
00244 *> \param[out] V
00245 *> \verbatim
00246 *>          V is COMPLEX array, dimension (LDV,P)
00247 *>          If JOBV = 'V', V contains the P-by-P unitary matrix V.
00248 *>          If JOBV = 'N', V is not referenced.
00249 *> \endverbatim
00250 *>
00251 *> \param[in] LDV
00252 *> \verbatim
00253 *>          LDV is INTEGER
00254 *>          The leading dimension of the array V. LDV >= max(1,P) if
00255 *>          JOBV = 'V'; LDV >= 1 otherwise.
00256 *> \endverbatim
00257 *>
00258 *> \param[out] Q
00259 *> \verbatim
00260 *>          Q is COMPLEX array, dimension (LDQ,N)
00261 *>          If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q.
00262 *>          If JOBQ = 'N', Q is not referenced.
00263 *> \endverbatim
00264 *>
00265 *> \param[in] LDQ
00266 *> \verbatim
00267 *>          LDQ is INTEGER
00268 *>          The leading dimension of the array Q. LDQ >= max(1,N) if
00269 *>          JOBQ = 'Q'; LDQ >= 1 otherwise.
00270 *> \endverbatim
00271 *>
00272 *> \param[out] WORK
00273 *> \verbatim
00274 *>          WORK is COMPLEX array, dimension (max(3*N,M,P)+N)
00275 *> \endverbatim
00276 *>
00277 *> \param[out] RWORK
00278 *> \verbatim
00279 *>          RWORK is REAL array, dimension (2*N)
00280 *> \endverbatim
00281 *>
00282 *> \param[out] IWORK
00283 *> \verbatim
00284 *>          IWORK is INTEGER array, dimension (N)
00285 *>          On exit, IWORK stores the sorting information. More
00286 *>          precisely, the following loop will sort ALPHA
00287 *>             for I = K+1, min(M,K+L)
00288 *>                 swap ALPHA(I) and ALPHA(IWORK(I))
00289 *>             endfor
00290 *>          such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N).
00291 *> \endverbatim
00292 *>
00293 *> \param[out] INFO
00294 *> \verbatim
00295 *>          INFO is INTEGER
00296 *>          = 0:  successful exit.
00297 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00298 *>          > 0:  if INFO = 1, the Jacobi-type procedure failed to
00299 *>                converge.  For further details, see subroutine CTGSJA.
00300 *> \endverbatim
00301 *
00302 *> \par Internal Parameters:
00303 *  =========================
00304 *>
00305 *> \verbatim
00306 *>  TOLA    REAL
00307 *>  TOLB    REAL
00308 *>          TOLA and TOLB are the thresholds to determine the effective
00309 *>          rank of (A**H,B**H)**H. Generally, they are set to
00310 *>                   TOLA = MAX(M,N)*norm(A)*MACHEPS,
00311 *>                   TOLB = MAX(P,N)*norm(B)*MACHEPS.
00312 *>          The size of TOLA and TOLB may affect the size of backward
00313 *>          errors of the decomposition.
00314 *> \endverbatim
00315 *
00316 *  Authors:
00317 *  ========
00318 *
00319 *> \author Univ. of Tennessee 
00320 *> \author Univ. of California Berkeley 
00321 *> \author Univ. of Colorado Denver 
00322 *> \author NAG Ltd. 
00323 *
00324 *> \date November 2011
00325 *
00326 *> \ingroup complexOTHERsing
00327 *
00328 *> \par Contributors:
00329 *  ==================
00330 *>
00331 *>     Ming Gu and Huan Ren, Computer Science Division, University of
00332 *>     California at Berkeley, USA
00333 *>
00334 *  =====================================================================
00335       SUBROUTINE CGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
00336      $                   LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK,
00337      $                   RWORK, IWORK, INFO )
00338 *
00339 *  -- LAPACK driver routine (version 3.4.0) --
00340 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00341 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00342 *     November 2011
00343 *
00344 *     .. Scalar Arguments ..
00345       CHARACTER          JOBQ, JOBU, JOBV
00346       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
00347 *     ..
00348 *     .. Array Arguments ..
00349       INTEGER            IWORK( * )
00350       REAL               ALPHA( * ), BETA( * ), RWORK( * )
00351       COMPLEX            A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
00352      $                   U( LDU, * ), V( LDV, * ), WORK( * )
00353 *     ..
00354 *
00355 *  =====================================================================
00356 *
00357 *     .. Local Scalars ..
00358       LOGICAL            WANTQ, WANTU, WANTV
00359       INTEGER            I, IBND, ISUB, J, NCYCLE
00360       REAL               ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
00361 *     ..
00362 *     .. External Functions ..
00363       LOGICAL            LSAME
00364       REAL               CLANGE, SLAMCH
00365       EXTERNAL           LSAME, CLANGE, SLAMCH
00366 *     ..
00367 *     .. External Subroutines ..
00368       EXTERNAL           CGGSVP, CTGSJA, SCOPY, XERBLA
00369 *     ..
00370 *     .. Intrinsic Functions ..
00371       INTRINSIC          MAX, MIN
00372 *     ..
00373 *     .. Executable Statements ..
00374 *
00375 *     Decode and test the input parameters
00376 *
00377       WANTU = LSAME( JOBU, 'U' )
00378       WANTV = LSAME( JOBV, 'V' )
00379       WANTQ = LSAME( JOBQ, 'Q' )
00380 *
00381       INFO = 0
00382       IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN
00383          INFO = -1
00384       ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN
00385          INFO = -2
00386       ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN
00387          INFO = -3
00388       ELSE IF( M.LT.0 ) THEN
00389          INFO = -4
00390       ELSE IF( N.LT.0 ) THEN
00391          INFO = -5
00392       ELSE IF( P.LT.0 ) THEN
00393          INFO = -6
00394       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00395          INFO = -10
00396       ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
00397          INFO = -12
00398       ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
00399          INFO = -16
00400       ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
00401          INFO = -18
00402       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
00403          INFO = -20
00404       END IF
00405       IF( INFO.NE.0 ) THEN
00406          CALL XERBLA( 'CGGSVD', -INFO )
00407          RETURN
00408       END IF
00409 *
00410 *     Compute the Frobenius norm of matrices A and B
00411 *
00412       ANORM = CLANGE( '1', M, N, A, LDA, RWORK )
00413       BNORM = CLANGE( '1', P, N, B, LDB, RWORK )
00414 *
00415 *     Get machine precision and set up threshold for determining
00416 *     the effective numerical rank of the matrices A and B.
00417 *
00418       ULP = SLAMCH( 'Precision' )
00419       UNFL = SLAMCH( 'Safe Minimum' )
00420       TOLA = MAX( M, N )*MAX( ANORM, UNFL )*ULP
00421       TOLB = MAX( P, N )*MAX( BNORM, UNFL )*ULP
00422 *
00423       CALL CGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA,
00424      $             TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK,
00425      $             WORK, WORK( N+1 ), INFO )
00426 *
00427 *     Compute the GSVD of two upper "triangular" matrices
00428 *
00429       CALL CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB,
00430      $             TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ,
00431      $             WORK, NCYCLE, INFO )
00432 *
00433 *     Sort the singular values and store the pivot indices in IWORK
00434 *     Copy ALPHA to RWORK, then sort ALPHA in RWORK
00435 *
00436       CALL SCOPY( N, ALPHA, 1, RWORK, 1 )
00437       IBND = MIN( L, M-K )
00438       DO 20 I = 1, IBND
00439 *
00440 *        Scan for largest ALPHA(K+I)
00441 *
00442          ISUB = I
00443          SMAX = RWORK( K+I )
00444          DO 10 J = I + 1, IBND
00445             TEMP = RWORK( K+J )
00446             IF( TEMP.GT.SMAX ) THEN
00447                ISUB = J
00448                SMAX = TEMP
00449             END IF
00450    10    CONTINUE
00451          IF( ISUB.NE.I ) THEN
00452             RWORK( K+ISUB ) = RWORK( K+I )
00453             RWORK( K+I ) = SMAX
00454             IWORK( K+I ) = K + ISUB
00455          ELSE
00456             IWORK( K+I ) = K + I
00457          END IF
00458    20 CONTINUE
00459 *
00460       RETURN
00461 *
00462 *     End of CGGSVD
00463 *
00464       END
 All Files Functions