![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b STGSJA 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download STGSJA + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgsja.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgsja.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgsja.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, 00022 * LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, 00023 * Q, LDQ, WORK, NCYCLE, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER JOBQ, JOBU, JOBV 00027 * INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, 00028 * $ NCYCLE, P 00029 * REAL TOLA, TOLB 00030 * .. 00031 * .. Array Arguments .. 00032 * REAL A( LDA, * ), ALPHA( * ), B( LDB, * ), 00033 * $ BETA( * ), Q( LDQ, * ), U( LDU, * ), 00034 * $ V( LDV, * ), WORK( * ) 00035 * .. 00036 * 00037 * 00038 *> \par Purpose: 00039 * ============= 00040 *> 00041 *> \verbatim 00042 *> 00043 *> STGSJA computes the generalized singular value decomposition (GSVD) 00044 *> of two real upper triangular (or trapezoidal) matrices A and B. 00045 *> 00046 *> On entry, it is assumed that matrices A and B have the following 00047 *> forms, which may be obtained by the preprocessing subroutine SGGSVP 00048 *> from a general M-by-N matrix A and P-by-N matrix B: 00049 *> 00050 *> N-K-L K L 00051 *> A = K ( 0 A12 A13 ) if M-K-L >= 0; 00052 *> L ( 0 0 A23 ) 00053 *> M-K-L ( 0 0 0 ) 00054 *> 00055 *> N-K-L K L 00056 *> A = K ( 0 A12 A13 ) if M-K-L < 0; 00057 *> M-K ( 0 0 A23 ) 00058 *> 00059 *> N-K-L K L 00060 *> B = L ( 0 0 B13 ) 00061 *> P-L ( 0 0 0 ) 00062 *> 00063 *> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular 00064 *> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, 00065 *> otherwise A23 is (M-K)-by-L upper trapezoidal. 00066 *> 00067 *> On exit, 00068 *> 00069 *> U**T *A*Q = D1*( 0 R ), V**T *B*Q = D2*( 0 R ), 00070 *> 00071 *> where U, V and Q are orthogonal matrices. 00072 *> R is a nonsingular upper triangular matrix, and D1 and D2 are 00073 *> ``diagonal'' matrices, which are of the following structures: 00074 *> 00075 *> If M-K-L >= 0, 00076 *> 00077 *> K L 00078 *> D1 = K ( I 0 ) 00079 *> L ( 0 C ) 00080 *> M-K-L ( 0 0 ) 00081 *> 00082 *> K L 00083 *> D2 = L ( 0 S ) 00084 *> P-L ( 0 0 ) 00085 *> 00086 *> N-K-L K L 00087 *> ( 0 R ) = K ( 0 R11 R12 ) K 00088 *> L ( 0 0 R22 ) L 00089 *> 00090 *> where 00091 *> 00092 *> C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), 00093 *> S = diag( BETA(K+1), ... , BETA(K+L) ), 00094 *> C**2 + S**2 = I. 00095 *> 00096 *> R is stored in A(1:K+L,N-K-L+1:N) on exit. 00097 *> 00098 *> If M-K-L < 0, 00099 *> 00100 *> K M-K K+L-M 00101 *> D1 = K ( I 0 0 ) 00102 *> M-K ( 0 C 0 ) 00103 *> 00104 *> K M-K K+L-M 00105 *> D2 = M-K ( 0 S 0 ) 00106 *> K+L-M ( 0 0 I ) 00107 *> P-L ( 0 0 0 ) 00108 *> 00109 *> N-K-L K M-K K+L-M 00110 *> ( 0 R ) = K ( 0 R11 R12 R13 ) 00111 *> M-K ( 0 0 R22 R23 ) 00112 *> K+L-M ( 0 0 0 R33 ) 00113 *> 00114 *> where 00115 *> C = diag( ALPHA(K+1), ... , ALPHA(M) ), 00116 *> S = diag( BETA(K+1), ... , BETA(M) ), 00117 *> C**2 + S**2 = I. 00118 *> 00119 *> R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored 00120 *> ( 0 R22 R23 ) 00121 *> in B(M-K+1:L,N+M-K-L+1:N) on exit. 00122 *> 00123 *> The computation of the orthogonal transformation matrices U, V or Q 00124 *> is optional. These matrices may either be formed explicitly, or they 00125 *> may be postmultiplied into input matrices U1, V1, or Q1. 00126 *> \endverbatim 00127 * 00128 * Arguments: 00129 * ========== 00130 * 00131 *> \param[in] JOBU 00132 *> \verbatim 00133 *> JOBU is CHARACTER*1 00134 *> = 'U': U must contain an orthogonal matrix U1 on entry, and 00135 *> the product U1*U is returned; 00136 *> = 'I': U is initialized to the unit matrix, and the 00137 *> orthogonal matrix U is returned; 00138 *> = 'N': U is not computed. 00139 *> \endverbatim 00140 *> 00141 *> \param[in] JOBV 00142 *> \verbatim 00143 *> JOBV is CHARACTER*1 00144 *> = 'V': V must contain an orthogonal matrix V1 on entry, and 00145 *> the product V1*V is returned; 00146 *> = 'I': V is initialized to the unit matrix, and the 00147 *> orthogonal matrix V is returned; 00148 *> = 'N': V is not computed. 00149 *> \endverbatim 00150 *> 00151 *> \param[in] JOBQ 00152 *> \verbatim 00153 *> JOBQ is CHARACTER*1 00154 *> = 'Q': Q must contain an orthogonal matrix Q1 on entry, and 00155 *> the product Q1*Q is returned; 00156 *> = 'I': Q is initialized to the unit matrix, and the 00157 *> orthogonal matrix Q is returned; 00158 *> = 'N': Q is not computed. 00159 *> \endverbatim 00160 *> 00161 *> \param[in] M 00162 *> \verbatim 00163 *> M is INTEGER 00164 *> The number of rows of the matrix A. M >= 0. 00165 *> \endverbatim 00166 *> 00167 *> \param[in] P 00168 *> \verbatim 00169 *> P is INTEGER 00170 *> The number of rows of the matrix B. P >= 0. 00171 *> \endverbatim 00172 *> 00173 *> \param[in] N 00174 *> \verbatim 00175 *> N is INTEGER 00176 *> The number of columns of the matrices A and B. N >= 0. 00177 *> \endverbatim 00178 *> 00179 *> \param[in] K 00180 *> \verbatim 00181 *> K is INTEGER 00182 *> \endverbatim 00183 *> 00184 *> \param[in] L 00185 *> \verbatim 00186 *> L is INTEGER 00187 *> 00188 *> K and L specify the subblocks in the input matrices A and B: 00189 *> A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,N-L+1:N) 00190 *> of A and B, whose GSVD is going to be computed by STGSJA. 00191 *> See Further Details. 00192 *> \endverbatim 00193 *> 00194 *> \param[in,out] A 00195 *> \verbatim 00196 *> A is REAL array, dimension (LDA,N) 00197 *> On entry, the M-by-N matrix A. 00198 *> On exit, A(N-K+1:N,1:MIN(K+L,M) ) contains the triangular 00199 *> matrix R or part of R. See Purpose for details. 00200 *> \endverbatim 00201 *> 00202 *> \param[in] LDA 00203 *> \verbatim 00204 *> LDA is INTEGER 00205 *> The leading dimension of the array A. LDA >= max(1,M). 00206 *> \endverbatim 00207 *> 00208 *> \param[in,out] B 00209 *> \verbatim 00210 *> B is REAL array, dimension (LDB,N) 00211 *> On entry, the P-by-N matrix B. 00212 *> On exit, if necessary, B(M-K+1:L,N+M-K-L+1:N) contains 00213 *> a part of R. See Purpose for details. 00214 *> \endverbatim 00215 *> 00216 *> \param[in] LDB 00217 *> \verbatim 00218 *> LDB is INTEGER 00219 *> The leading dimension of the array B. LDB >= max(1,P). 00220 *> \endverbatim 00221 *> 00222 *> \param[in] TOLA 00223 *> \verbatim 00224 *> TOLA is REAL 00225 *> \endverbatim 00226 *> 00227 *> \param[in] TOLB 00228 *> \verbatim 00229 *> TOLB is REAL 00230 *> 00231 *> TOLA and TOLB are the convergence criteria for the Jacobi- 00232 *> Kogbetliantz iteration procedure. Generally, they are the 00233 *> same as used in the preprocessing step, say 00234 *> TOLA = max(M,N)*norm(A)*MACHEPS, 00235 *> TOLB = max(P,N)*norm(B)*MACHEPS. 00236 *> \endverbatim 00237 *> 00238 *> \param[out] ALPHA 00239 *> \verbatim 00240 *> ALPHA is REAL array, dimension (N) 00241 *> \endverbatim 00242 *> 00243 *> \param[out] BETA 00244 *> \verbatim 00245 *> BETA is REAL array, dimension (N) 00246 *> 00247 *> On exit, ALPHA and BETA contain the generalized singular 00248 *> value pairs of A and B; 00249 *> ALPHA(1:K) = 1, 00250 *> BETA(1:K) = 0, 00251 *> and if M-K-L >= 0, 00252 *> ALPHA(K+1:K+L) = diag(C), 00253 *> BETA(K+1:K+L) = diag(S), 00254 *> or if M-K-L < 0, 00255 *> ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0 00256 *> BETA(K+1:M) = S, BETA(M+1:K+L) = 1. 00257 *> Furthermore, if K+L < N, 00258 *> ALPHA(K+L+1:N) = 0 and 00259 *> BETA(K+L+1:N) = 0. 00260 *> \endverbatim 00261 *> 00262 *> \param[in,out] U 00263 *> \verbatim 00264 *> U is REAL array, dimension (LDU,M) 00265 *> On entry, if JOBU = 'U', U must contain a matrix U1 (usually 00266 *> the orthogonal matrix returned by SGGSVP). 00267 *> On exit, 00268 *> if JOBU = 'I', U contains the orthogonal matrix U; 00269 *> if JOBU = 'U', U contains the product U1*U. 00270 *> If JOBU = 'N', U is not referenced. 00271 *> \endverbatim 00272 *> 00273 *> \param[in] LDU 00274 *> \verbatim 00275 *> LDU is INTEGER 00276 *> The leading dimension of the array U. LDU >= max(1,M) if 00277 *> JOBU = 'U'; LDU >= 1 otherwise. 00278 *> \endverbatim 00279 *> 00280 *> \param[in,out] V 00281 *> \verbatim 00282 *> V is REAL array, dimension (LDV,P) 00283 *> On entry, if JOBV = 'V', V must contain a matrix V1 (usually 00284 *> the orthogonal matrix returned by SGGSVP). 00285 *> On exit, 00286 *> if JOBV = 'I', V contains the orthogonal matrix V; 00287 *> if JOBV = 'V', V contains the product V1*V. 00288 *> If JOBV = 'N', V is not referenced. 00289 *> \endverbatim 00290 *> 00291 *> \param[in] LDV 00292 *> \verbatim 00293 *> LDV is INTEGER 00294 *> The leading dimension of the array V. LDV >= max(1,P) if 00295 *> JOBV = 'V'; LDV >= 1 otherwise. 00296 *> \endverbatim 00297 *> 00298 *> \param[in,out] Q 00299 *> \verbatim 00300 *> Q is REAL array, dimension (LDQ,N) 00301 *> On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually 00302 *> the orthogonal matrix returned by SGGSVP). 00303 *> On exit, 00304 *> if JOBQ = 'I', Q contains the orthogonal matrix Q; 00305 *> if JOBQ = 'Q', Q contains the product Q1*Q. 00306 *> If JOBQ = 'N', Q is not referenced. 00307 *> \endverbatim 00308 *> 00309 *> \param[in] LDQ 00310 *> \verbatim 00311 *> LDQ is INTEGER 00312 *> The leading dimension of the array Q. LDQ >= max(1,N) if 00313 *> JOBQ = 'Q'; LDQ >= 1 otherwise. 00314 *> \endverbatim 00315 *> 00316 *> \param[out] WORK 00317 *> \verbatim 00318 *> WORK is REAL array, dimension (2*N) 00319 *> \endverbatim 00320 *> 00321 *> \param[out] NCYCLE 00322 *> \verbatim 00323 *> NCYCLE is INTEGER 00324 *> The number of cycles required for convergence. 00325 *> \endverbatim 00326 *> 00327 *> \param[out] INFO 00328 *> \verbatim 00329 *> INFO is INTEGER 00330 *> = 0: successful exit 00331 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00332 *> = 1: the procedure does not converge after MAXIT cycles. 00333 *> \endverbatim 00334 *> 00335 *> \verbatim 00336 *> Internal Parameters 00337 *> =================== 00338 *> 00339 *> MAXIT INTEGER 00340 *> MAXIT specifies the total loops that the iterative procedure 00341 *> may take. If after MAXIT cycles, the routine fails to 00342 *> converge, we return INFO = 1. 00343 *> \endverbatim 00344 * 00345 * Authors: 00346 * ======== 00347 * 00348 *> \author Univ. of Tennessee 00349 *> \author Univ. of California Berkeley 00350 *> \author Univ. of Colorado Denver 00351 *> \author NAG Ltd. 00352 * 00353 *> \date November 2011 00354 * 00355 *> \ingroup realOTHERcomputational 00356 * 00357 *> \par Further Details: 00358 * ===================== 00359 *> 00360 *> \verbatim 00361 *> 00362 *> STGSJA essentially uses a variant of Kogbetliantz algorithm to reduce 00363 *> min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L 00364 *> matrix B13 to the form: 00365 *> 00366 *> U1**T *A13*Q1 = C1*R1; V1**T *B13*Q1 = S1*R1, 00367 *> 00368 *> where U1, V1 and Q1 are orthogonal matrix, and Z**T is the transpose 00369 *> of Z. C1 and S1 are diagonal matrices satisfying 00370 *> 00371 *> C1**2 + S1**2 = I, 00372 *> 00373 *> and R1 is an L-by-L nonsingular upper triangular matrix. 00374 *> \endverbatim 00375 *> 00376 * ===================================================================== 00377 SUBROUTINE STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, 00378 $ LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, 00379 $ Q, LDQ, WORK, NCYCLE, INFO ) 00380 * 00381 * -- LAPACK computational routine (version 3.4.0) -- 00382 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00383 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00384 * November 2011 00385 * 00386 * .. Scalar Arguments .. 00387 CHARACTER JOBQ, JOBU, JOBV 00388 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, 00389 $ NCYCLE, P 00390 REAL TOLA, TOLB 00391 * .. 00392 * .. Array Arguments .. 00393 REAL A( LDA, * ), ALPHA( * ), B( LDB, * ), 00394 $ BETA( * ), Q( LDQ, * ), U( LDU, * ), 00395 $ V( LDV, * ), WORK( * ) 00396 * .. 00397 * 00398 * ===================================================================== 00399 * 00400 * .. Parameters .. 00401 INTEGER MAXIT 00402 PARAMETER ( MAXIT = 40 ) 00403 REAL ZERO, ONE 00404 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00405 * .. 00406 * .. Local Scalars .. 00407 * 00408 LOGICAL INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV 00409 INTEGER I, J, KCYCLE 00410 REAL A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, ERROR, 00411 $ GAMMA, RWK, SNQ, SNU, SNV, SSMIN 00412 * .. 00413 * .. External Functions .. 00414 LOGICAL LSAME 00415 EXTERNAL LSAME 00416 * .. 00417 * .. External Subroutines .. 00418 EXTERNAL SCOPY, SLAGS2, SLAPLL, SLARTG, SLASET, SROT, 00419 $ SSCAL, XERBLA 00420 * .. 00421 * .. Intrinsic Functions .. 00422 INTRINSIC ABS, MAX, MIN 00423 * .. 00424 * .. Executable Statements .. 00425 * 00426 * Decode and test the input parameters 00427 * 00428 INITU = LSAME( JOBU, 'I' ) 00429 WANTU = INITU .OR. LSAME( JOBU, 'U' ) 00430 * 00431 INITV = LSAME( JOBV, 'I' ) 00432 WANTV = INITV .OR. LSAME( JOBV, 'V' ) 00433 * 00434 INITQ = LSAME( JOBQ, 'I' ) 00435 WANTQ = INITQ .OR. LSAME( JOBQ, 'Q' ) 00436 * 00437 INFO = 0 00438 IF( .NOT.( INITU .OR. WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN 00439 INFO = -1 00440 ELSE IF( .NOT.( INITV .OR. WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN 00441 INFO = -2 00442 ELSE IF( .NOT.( INITQ .OR. WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN 00443 INFO = -3 00444 ELSE IF( M.LT.0 ) THEN 00445 INFO = -4 00446 ELSE IF( P.LT.0 ) THEN 00447 INFO = -5 00448 ELSE IF( N.LT.0 ) THEN 00449 INFO = -6 00450 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00451 INFO = -10 00452 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 00453 INFO = -12 00454 ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN 00455 INFO = -18 00456 ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN 00457 INFO = -20 00458 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 00459 INFO = -22 00460 END IF 00461 IF( INFO.NE.0 ) THEN 00462 CALL XERBLA( 'STGSJA', -INFO ) 00463 RETURN 00464 END IF 00465 * 00466 * Initialize U, V and Q, if necessary 00467 * 00468 IF( INITU ) 00469 $ CALL SLASET( 'Full', M, M, ZERO, ONE, U, LDU ) 00470 IF( INITV ) 00471 $ CALL SLASET( 'Full', P, P, ZERO, ONE, V, LDV ) 00472 IF( INITQ ) 00473 $ CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) 00474 * 00475 * Loop until convergence 00476 * 00477 UPPER = .FALSE. 00478 DO 40 KCYCLE = 1, MAXIT 00479 * 00480 UPPER = .NOT.UPPER 00481 * 00482 DO 20 I = 1, L - 1 00483 DO 10 J = I + 1, L 00484 * 00485 A1 = ZERO 00486 A2 = ZERO 00487 A3 = ZERO 00488 IF( K+I.LE.M ) 00489 $ A1 = A( K+I, N-L+I ) 00490 IF( K+J.LE.M ) 00491 $ A3 = A( K+J, N-L+J ) 00492 * 00493 B1 = B( I, N-L+I ) 00494 B3 = B( J, N-L+J ) 00495 * 00496 IF( UPPER ) THEN 00497 IF( K+I.LE.M ) 00498 $ A2 = A( K+I, N-L+J ) 00499 B2 = B( I, N-L+J ) 00500 ELSE 00501 IF( K+J.LE.M ) 00502 $ A2 = A( K+J, N-L+I ) 00503 B2 = B( J, N-L+I ) 00504 END IF 00505 * 00506 CALL SLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, 00507 $ CSV, SNV, CSQ, SNQ ) 00508 * 00509 * Update (K+I)-th and (K+J)-th rows of matrix A: U**T *A 00510 * 00511 IF( K+J.LE.M ) 00512 $ CALL SROT( L, A( K+J, N-L+1 ), LDA, A( K+I, N-L+1 ), 00513 $ LDA, CSU, SNU ) 00514 * 00515 * Update I-th and J-th rows of matrix B: V**T *B 00516 * 00517 CALL SROT( L, B( J, N-L+1 ), LDB, B( I, N-L+1 ), LDB, 00518 $ CSV, SNV ) 00519 * 00520 * Update (N-L+I)-th and (N-L+J)-th columns of matrices 00521 * A and B: A*Q and B*Q 00522 * 00523 CALL SROT( MIN( K+L, M ), A( 1, N-L+J ), 1, 00524 $ A( 1, N-L+I ), 1, CSQ, SNQ ) 00525 * 00526 CALL SROT( L, B( 1, N-L+J ), 1, B( 1, N-L+I ), 1, CSQ, 00527 $ SNQ ) 00528 * 00529 IF( UPPER ) THEN 00530 IF( K+I.LE.M ) 00531 $ A( K+I, N-L+J ) = ZERO 00532 B( I, N-L+J ) = ZERO 00533 ELSE 00534 IF( K+J.LE.M ) 00535 $ A( K+J, N-L+I ) = ZERO 00536 B( J, N-L+I ) = ZERO 00537 END IF 00538 * 00539 * Update orthogonal matrices U, V, Q, if desired. 00540 * 00541 IF( WANTU .AND. K+J.LE.M ) 00542 $ CALL SROT( M, U( 1, K+J ), 1, U( 1, K+I ), 1, CSU, 00543 $ SNU ) 00544 * 00545 IF( WANTV ) 00546 $ CALL SROT( P, V( 1, J ), 1, V( 1, I ), 1, CSV, SNV ) 00547 * 00548 IF( WANTQ ) 00549 $ CALL SROT( N, Q( 1, N-L+J ), 1, Q( 1, N-L+I ), 1, CSQ, 00550 $ SNQ ) 00551 * 00552 10 CONTINUE 00553 20 CONTINUE 00554 * 00555 IF( .NOT.UPPER ) THEN 00556 * 00557 * The matrices A13 and B13 were lower triangular at the start 00558 * of the cycle, and are now upper triangular. 00559 * 00560 * Convergence test: test the parallelism of the corresponding 00561 * rows of A and B. 00562 * 00563 ERROR = ZERO 00564 DO 30 I = 1, MIN( L, M-K ) 00565 CALL SCOPY( L-I+1, A( K+I, N-L+I ), LDA, WORK, 1 ) 00566 CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, WORK( L+1 ), 1 ) 00567 CALL SLAPLL( L-I+1, WORK, 1, WORK( L+1 ), 1, SSMIN ) 00568 ERROR = MAX( ERROR, SSMIN ) 00569 30 CONTINUE 00570 * 00571 IF( ABS( ERROR ).LE.MIN( TOLA, TOLB ) ) 00572 $ GO TO 50 00573 END IF 00574 * 00575 * End of cycle loop 00576 * 00577 40 CONTINUE 00578 * 00579 * The algorithm has not converged after MAXIT cycles. 00580 * 00581 INFO = 1 00582 GO TO 100 00583 * 00584 50 CONTINUE 00585 * 00586 * If ERROR <= MIN(TOLA,TOLB), then the algorithm has converged. 00587 * Compute the generalized singular value pairs (ALPHA, BETA), and 00588 * set the triangular matrix R to array A. 00589 * 00590 DO 60 I = 1, K 00591 ALPHA( I ) = ONE 00592 BETA( I ) = ZERO 00593 60 CONTINUE 00594 * 00595 DO 70 I = 1, MIN( L, M-K ) 00596 * 00597 A1 = A( K+I, N-L+I ) 00598 B1 = B( I, N-L+I ) 00599 * 00600 IF( A1.NE.ZERO ) THEN 00601 GAMMA = B1 / A1 00602 * 00603 * change sign if necessary 00604 * 00605 IF( GAMMA.LT.ZERO ) THEN 00606 CALL SSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB ) 00607 IF( WANTV ) 00608 $ CALL SSCAL( P, -ONE, V( 1, I ), 1 ) 00609 END IF 00610 * 00611 CALL SLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ), 00612 $ RWK ) 00613 * 00614 IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN 00615 CALL SSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ), 00616 $ LDA ) 00617 ELSE 00618 CALL SSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ), 00619 $ LDB ) 00620 CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ), 00621 $ LDA ) 00622 END IF 00623 * 00624 ELSE 00625 * 00626 ALPHA( K+I ) = ZERO 00627 BETA( K+I ) = ONE 00628 CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ), 00629 $ LDA ) 00630 * 00631 END IF 00632 * 00633 70 CONTINUE 00634 * 00635 * Post-assignment 00636 * 00637 DO 80 I = M + 1, K + L 00638 ALPHA( I ) = ZERO 00639 BETA( I ) = ONE 00640 80 CONTINUE 00641 * 00642 IF( K+L.LT.N ) THEN 00643 DO 90 I = K + L + 1, N 00644 ALPHA( I ) = ZERO 00645 BETA( I ) = ZERO 00646 90 CONTINUE 00647 END IF 00648 * 00649 100 CONTINUE 00650 NCYCLE = KCYCLE 00651 RETURN 00652 * 00653 * End of STGSJA 00654 * 00655 END