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