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