LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sggsvd.f
Go to the documentation of this file.
00001 *> \brief <b> SGGSVD 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 SGGSVD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sggsvd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sggsvd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sggsvd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
00022 *                          LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK,
00023 *                          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               A( LDA, * ), ALPHA( * ), B( LDB, * ),
00032 *      $                   BETA( * ), Q( LDQ, * ), U( LDU, * ),
00033 *      $                   V( LDV, * ), WORK( * )
00034 *       ..
00035 *  
00036 *
00037 *> \par Purpose:
00038 *  =============
00039 *>
00040 *> \verbatim
00041 *>
00042 *> SGGSVD computes the generalized singular value decomposition (GSVD)
00043 *> of an M-by-N real matrix A and P-by-N real matrix B:
00044 *>
00045 *>       U**T*A*Q = D1*( 0 R ),    V**T*B*Q = D2*( 0 R )
00046 *>
00047 *> where U, V and Q are orthogonal matrices.
00048 *> Let K+L = the effective numerical rank of the matrix (A**T,B**T)**T,
00049 *> then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and
00050 *> D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the
00051 *> 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 orthogonal
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**T.
00108 *> If ( A**T,B**T)**T  has orthonormal columns, then the GSVD of A and B is
00109 *> also equal to the CS decomposition of A and B. Furthermore, the GSVD
00110 *> can be used to derive the solution of the eigenvalue problem:
00111 *>                      A**T*A x = lambda* B**T*B x.
00112 *> In some literature, the GSVD of A and B is presented in the form
00113 *>                  U**T*A*X = ( 0 D1 ),   V**T*B*X = ( 0 D2 )
00114 *> where U and V are orthogonal and X is nonsingular, 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':  Orthogonal 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':  Orthogonal 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':  Orthogonal 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**T,B**T)**T.
00176 *> \endverbatim
00177 *>
00178 *> \param[in,out] A
00179 *> \verbatim
00180 *>          A is REAL 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 REAL array, dimension (LDB,N)
00195 *>          On entry, the P-by-N matrix B.
00196 *>          On exit, B contains the triangular matrix R if M-K-L < 0.
00197 *>          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 REAL array, dimension (LDU,M)
00233 *>          If JOBU = 'U', U contains the M-by-M orthogonal 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 REAL array, dimension (LDV,P)
00247 *>          If JOBV = 'V', V contains the P-by-P orthogonal 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 REAL array, dimension (LDQ,N)
00261 *>          If JOBQ = 'Q', Q contains the N-by-N orthogonal 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 REAL array,
00275 *>                      dimension (max(3*N,M,P)+N)
00276 *> \endverbatim
00277 *>
00278 *> \param[out] IWORK
00279 *> \verbatim
00280 *>          IWORK is INTEGER array, dimension (N)
00281 *>          On exit, IWORK stores the sorting information. More
00282 *>          precisely, the following loop will sort ALPHA
00283 *>             for I = K+1, min(M,K+L)
00284 *>                 swap ALPHA(I) and ALPHA(IWORK(I))
00285 *>             endfor
00286 *>          such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N).
00287 *> \endverbatim
00288 *>
00289 *> \param[out] INFO
00290 *> \verbatim
00291 *>          INFO is INTEGER
00292 *>          = 0:  successful exit
00293 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00294 *>          > 0:  if INFO = 1, the Jacobi-type procedure failed to
00295 *>                converge.  For further details, see subroutine STGSJA.
00296 *> \endverbatim
00297 *
00298 *> \par Internal Parameters:
00299 *  =========================
00300 *>
00301 *> \verbatim
00302 *>  TOLA    REAL
00303 *>  TOLB    REAL
00304 *>          TOLA and TOLB are the thresholds to determine the effective
00305 *>          rank of (A**T,B**T)**T. Generally, they are set to
00306 *>                   TOLA = MAX(M,N)*norm(A)*MACHEPS,
00307 *>                   TOLB = MAX(P,N)*norm(B)*MACHEPS.
00308 *>          The size of TOLA and TOLB may affect the size of backward
00309 *>          errors of the decomposition.
00310 *> \endverbatim
00311 *
00312 *  Authors:
00313 *  ========
00314 *
00315 *> \author Univ. of Tennessee 
00316 *> \author Univ. of California Berkeley 
00317 *> \author Univ. of Colorado Denver 
00318 *> \author NAG Ltd. 
00319 *
00320 *> \date November 2011
00321 *
00322 *> \ingroup realOTHERsing
00323 *
00324 *> \par Contributors:
00325 *  ==================
00326 *>
00327 *>     Ming Gu and Huan Ren, Computer Science Division, University of
00328 *>     California at Berkeley, USA
00329 *>
00330 *  =====================================================================
00331       SUBROUTINE SGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
00332      $                   LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK,
00333      $                   IWORK, INFO )
00334 *
00335 *  -- LAPACK driver routine (version 3.4.0) --
00336 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00337 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00338 *     November 2011
00339 *
00340 *     .. Scalar Arguments ..
00341       CHARACTER          JOBQ, JOBU, JOBV
00342       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
00343 *     ..
00344 *     .. Array Arguments ..
00345       INTEGER            IWORK( * )
00346       REAL               A( LDA, * ), ALPHA( * ), B( LDB, * ),
00347      $                   BETA( * ), Q( LDQ, * ), U( LDU, * ),
00348      $                   V( LDV, * ), WORK( * )
00349 *     ..
00350 *
00351 *  =====================================================================
00352 *
00353 *     .. Local Scalars ..
00354       LOGICAL            WANTQ, WANTU, WANTV
00355       INTEGER            I, IBND, ISUB, J, NCYCLE
00356       REAL               ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
00357 *     ..
00358 *     .. External Functions ..
00359       LOGICAL            LSAME
00360       REAL               SLAMCH, SLANGE
00361       EXTERNAL           LSAME, SLAMCH, SLANGE
00362 *     ..
00363 *     .. External Subroutines ..
00364       EXTERNAL           SCOPY, SGGSVP, STGSJA, XERBLA
00365 *     ..
00366 *     .. Intrinsic Functions ..
00367       INTRINSIC          MAX, MIN
00368 *     ..
00369 *     .. Executable Statements ..
00370 *
00371 *     Test the input parameters
00372 *
00373       WANTU = LSAME( JOBU, 'U' )
00374       WANTV = LSAME( JOBV, 'V' )
00375       WANTQ = LSAME( JOBQ, 'Q' )
00376 *
00377       INFO = 0
00378       IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN
00379          INFO = -1
00380       ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN
00381          INFO = -2
00382       ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN
00383          INFO = -3
00384       ELSE IF( M.LT.0 ) THEN
00385          INFO = -4
00386       ELSE IF( N.LT.0 ) THEN
00387          INFO = -5
00388       ELSE IF( P.LT.0 ) THEN
00389          INFO = -6
00390       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00391          INFO = -10
00392       ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
00393          INFO = -12
00394       ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
00395          INFO = -16
00396       ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
00397          INFO = -18
00398       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
00399          INFO = -20
00400       END IF
00401       IF( INFO.NE.0 ) THEN
00402          CALL XERBLA( 'SGGSVD', -INFO )
00403          RETURN
00404       END IF
00405 *
00406 *     Compute the Frobenius norm of matrices A and B
00407 *
00408       ANORM = SLANGE( '1', M, N, A, LDA, WORK )
00409       BNORM = SLANGE( '1', P, N, B, LDB, WORK )
00410 *
00411 *     Get machine precision and set up threshold for determining
00412 *     the effective numerical rank of the matrices A and B.
00413 *
00414       ULP = SLAMCH( 'Precision' )
00415       UNFL = SLAMCH( 'Safe Minimum' )
00416       TOLA = MAX( M, N )*MAX( ANORM, UNFL )*ULP
00417       TOLB = MAX( P, N )*MAX( BNORM, UNFL )*ULP
00418 *
00419 *     Preprocessing
00420 *
00421       CALL SGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA,
00422      $             TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, WORK,
00423      $             WORK( N+1 ), INFO )
00424 *
00425 *     Compute the GSVD of two upper "triangular" matrices
00426 *
00427       CALL STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB,
00428      $             TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ,
00429      $             WORK, NCYCLE, INFO )
00430 *
00431 *     Sort the singular values and store the pivot indices in IWORK
00432 *     Copy ALPHA to WORK, then sort ALPHA in WORK
00433 *
00434       CALL SCOPY( N, ALPHA, 1, WORK, 1 )
00435       IBND = MIN( L, M-K )
00436       DO 20 I = 1, IBND
00437 *
00438 *        Scan for largest ALPHA(K+I)
00439 *
00440          ISUB = I
00441          SMAX = WORK( K+I )
00442          DO 10 J = I + 1, IBND
00443             TEMP = WORK( K+J )
00444             IF( TEMP.GT.SMAX ) THEN
00445                ISUB = J
00446                SMAX = TEMP
00447             END IF
00448    10    CONTINUE
00449          IF( ISUB.NE.I ) THEN
00450             WORK( K+ISUB ) = WORK( K+I )
00451             WORK( K+I ) = SMAX
00452             IWORK( K+I ) = K + ISUB
00453          ELSE
00454             IWORK( K+I ) = K + I
00455          END IF
00456    20 CONTINUE
00457 *
00458       RETURN
00459 *
00460 *     End of SGGSVD
00461 *
00462       END
 All Files Functions