![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SGGSVP 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SGGSVP + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sggsvp.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sggsvp.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sggsvp.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, 00022 * TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, 00023 * IWORK, TAU, WORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER JOBQ, JOBU, JOBV 00027 * INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P 00028 * REAL TOLA, TOLB 00029 * .. 00030 * .. Array Arguments .. 00031 * INTEGER IWORK( * ) 00032 * REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00033 * $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) 00034 * .. 00035 * 00036 * 00037 *> \par Purpose: 00038 * ============= 00039 *> 00040 *> \verbatim 00041 *> 00042 *> SGGSVP computes orthogonal matrices U, V and Q such that 00043 *> 00044 *> N-K-L K L 00045 *> U**T*A*Q = K ( 0 A12 A13 ) if M-K-L >= 0; 00046 *> L ( 0 0 A23 ) 00047 *> M-K-L ( 0 0 0 ) 00048 *> 00049 *> N-K-L K L 00050 *> = K ( 0 A12 A13 ) if M-K-L < 0; 00051 *> M-K ( 0 0 A23 ) 00052 *> 00053 *> N-K-L K L 00054 *> V**T*B*Q = L ( 0 0 B13 ) 00055 *> P-L ( 0 0 0 ) 00056 *> 00057 *> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular 00058 *> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, 00059 *> otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective 00060 *> numerical rank of the (M+P)-by-N matrix (A**T,B**T)**T. 00061 *> 00062 *> This decomposition is the preprocessing step for computing the 00063 *> Generalized Singular Value Decomposition (GSVD), see subroutine 00064 *> SGGSVD. 00065 *> \endverbatim 00066 * 00067 * Arguments: 00068 * ========== 00069 * 00070 *> \param[in] JOBU 00071 *> \verbatim 00072 *> JOBU is CHARACTER*1 00073 *> = 'U': Orthogonal matrix U is computed; 00074 *> = 'N': U is not computed. 00075 *> \endverbatim 00076 *> 00077 *> \param[in] JOBV 00078 *> \verbatim 00079 *> JOBV is CHARACTER*1 00080 *> = 'V': Orthogonal matrix V is computed; 00081 *> = 'N': V is not computed. 00082 *> \endverbatim 00083 *> 00084 *> \param[in] JOBQ 00085 *> \verbatim 00086 *> JOBQ is CHARACTER*1 00087 *> = 'Q': Orthogonal matrix Q is computed; 00088 *> = 'N': Q is not computed. 00089 *> \endverbatim 00090 *> 00091 *> \param[in] M 00092 *> \verbatim 00093 *> M is INTEGER 00094 *> The number of rows of the matrix A. M >= 0. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] P 00098 *> \verbatim 00099 *> P is INTEGER 00100 *> The number of rows of the matrix B. P >= 0. 00101 *> \endverbatim 00102 *> 00103 *> \param[in] N 00104 *> \verbatim 00105 *> N is INTEGER 00106 *> The number of columns of the matrices A and B. N >= 0. 00107 *> \endverbatim 00108 *> 00109 *> \param[in,out] A 00110 *> \verbatim 00111 *> A is REAL array, dimension (LDA,N) 00112 *> On entry, the M-by-N matrix A. 00113 *> On exit, A contains the triangular (or trapezoidal) matrix 00114 *> described in the Purpose section. 00115 *> \endverbatim 00116 *> 00117 *> \param[in] LDA 00118 *> \verbatim 00119 *> LDA is INTEGER 00120 *> The leading dimension of the array A. LDA >= max(1,M). 00121 *> \endverbatim 00122 *> 00123 *> \param[in,out] B 00124 *> \verbatim 00125 *> B is REAL array, dimension (LDB,N) 00126 *> On entry, the P-by-N matrix B. 00127 *> On exit, B contains the triangular matrix described in 00128 *> the Purpose section. 00129 *> \endverbatim 00130 *> 00131 *> \param[in] LDB 00132 *> \verbatim 00133 *> LDB is INTEGER 00134 *> The leading dimension of the array B. LDB >= max(1,P). 00135 *> \endverbatim 00136 *> 00137 *> \param[in] TOLA 00138 *> \verbatim 00139 *> TOLA is REAL 00140 *> \endverbatim 00141 *> 00142 *> \param[in] TOLB 00143 *> \verbatim 00144 *> TOLB is REAL 00145 *> 00146 *> TOLA and TOLB are the thresholds to determine the effective 00147 *> numerical rank of matrix B and a subblock of A. Generally, 00148 *> they are set to 00149 *> TOLA = MAX(M,N)*norm(A)*MACHEPS, 00150 *> TOLB = MAX(P,N)*norm(B)*MACHEPS. 00151 *> The size of TOLA and TOLB may affect the size of backward 00152 *> errors of the decomposition. 00153 *> \endverbatim 00154 *> 00155 *> \param[out] K 00156 *> \verbatim 00157 *> K is INTEGER 00158 *> \endverbatim 00159 *> 00160 *> \param[out] L 00161 *> \verbatim 00162 *> L is INTEGER 00163 *> 00164 *> On exit, K and L specify the dimension of the subblocks 00165 *> described in Purpose section. 00166 *> K + L = effective numerical rank of (A**T,B**T)**T. 00167 *> \endverbatim 00168 *> 00169 *> \param[out] U 00170 *> \verbatim 00171 *> U is REAL array, dimension (LDU,M) 00172 *> If JOBU = 'U', U contains the orthogonal matrix U. 00173 *> If JOBU = 'N', U is not referenced. 00174 *> \endverbatim 00175 *> 00176 *> \param[in] LDU 00177 *> \verbatim 00178 *> LDU is INTEGER 00179 *> The leading dimension of the array U. LDU >= max(1,M) if 00180 *> JOBU = 'U'; LDU >= 1 otherwise. 00181 *> \endverbatim 00182 *> 00183 *> \param[out] V 00184 *> \verbatim 00185 *> V is REAL array, dimension (LDV,P) 00186 *> If JOBV = 'V', V contains the orthogonal matrix V. 00187 *> If JOBV = 'N', V is not referenced. 00188 *> \endverbatim 00189 *> 00190 *> \param[in] LDV 00191 *> \verbatim 00192 *> LDV is INTEGER 00193 *> The leading dimension of the array V. LDV >= max(1,P) if 00194 *> JOBV = 'V'; LDV >= 1 otherwise. 00195 *> \endverbatim 00196 *> 00197 *> \param[out] Q 00198 *> \verbatim 00199 *> Q is REAL array, dimension (LDQ,N) 00200 *> If JOBQ = 'Q', Q contains the orthogonal matrix Q. 00201 *> If JOBQ = 'N', Q is not referenced. 00202 *> \endverbatim 00203 *> 00204 *> \param[in] LDQ 00205 *> \verbatim 00206 *> LDQ is INTEGER 00207 *> The leading dimension of the array Q. LDQ >= max(1,N) if 00208 *> JOBQ = 'Q'; LDQ >= 1 otherwise. 00209 *> \endverbatim 00210 *> 00211 *> \param[out] IWORK 00212 *> \verbatim 00213 *> IWORK is INTEGER array, dimension (N) 00214 *> \endverbatim 00215 *> 00216 *> \param[out] TAU 00217 *> \verbatim 00218 *> TAU is REAL array, dimension (N) 00219 *> \endverbatim 00220 *> 00221 *> \param[out] WORK 00222 *> \verbatim 00223 *> WORK is REAL array, dimension (max(3*N,M,P)) 00224 *> \endverbatim 00225 *> 00226 *> \param[out] INFO 00227 *> \verbatim 00228 *> INFO is INTEGER 00229 *> = 0: successful exit 00230 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00231 *> \endverbatim 00232 * 00233 * Authors: 00234 * ======== 00235 * 00236 *> \author Univ. of Tennessee 00237 *> \author Univ. of California Berkeley 00238 *> \author Univ. of Colorado Denver 00239 *> \author NAG Ltd. 00240 * 00241 *> \date November 2011 00242 * 00243 *> \ingroup realOTHERcomputational 00244 * 00245 *> \par Further Details: 00246 * ===================== 00247 *> 00248 *> The subroutine uses LAPACK subroutine SGEQPF for the QR factorization 00249 *> with column pivoting to detect the effective numerical rank of the 00250 *> a matrix. It may be replaced by a better rank determination strategy. 00251 *> 00252 * ===================================================================== 00253 SUBROUTINE SGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, 00254 $ TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, 00255 $ IWORK, TAU, WORK, INFO ) 00256 * 00257 * -- LAPACK computational routine (version 3.4.0) -- 00258 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00259 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00260 * November 2011 00261 * 00262 * .. Scalar Arguments .. 00263 CHARACTER JOBQ, JOBU, JOBV 00264 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P 00265 REAL TOLA, TOLB 00266 * .. 00267 * .. Array Arguments .. 00268 INTEGER IWORK( * ) 00269 REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00270 $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) 00271 * .. 00272 * 00273 * ===================================================================== 00274 * 00275 * .. Parameters .. 00276 REAL ZERO, ONE 00277 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00278 * .. 00279 * .. Local Scalars .. 00280 LOGICAL FORWRD, WANTQ, WANTU, WANTV 00281 INTEGER I, J 00282 * .. 00283 * .. External Functions .. 00284 LOGICAL LSAME 00285 EXTERNAL LSAME 00286 * .. 00287 * .. External Subroutines .. 00288 EXTERNAL SGEQPF, SGEQR2, SGERQ2, SLACPY, SLAPMT, SLASET, 00289 $ SORG2R, SORM2R, SORMR2, XERBLA 00290 * .. 00291 * .. Intrinsic Functions .. 00292 INTRINSIC ABS, MAX, MIN 00293 * .. 00294 * .. Executable Statements .. 00295 * 00296 * Test the input parameters 00297 * 00298 WANTU = LSAME( JOBU, 'U' ) 00299 WANTV = LSAME( JOBV, 'V' ) 00300 WANTQ = LSAME( JOBQ, 'Q' ) 00301 FORWRD = .TRUE. 00302 * 00303 INFO = 0 00304 IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN 00305 INFO = -1 00306 ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN 00307 INFO = -2 00308 ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN 00309 INFO = -3 00310 ELSE IF( M.LT.0 ) THEN 00311 INFO = -4 00312 ELSE IF( P.LT.0 ) THEN 00313 INFO = -5 00314 ELSE IF( N.LT.0 ) THEN 00315 INFO = -6 00316 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00317 INFO = -8 00318 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 00319 INFO = -10 00320 ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN 00321 INFO = -16 00322 ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN 00323 INFO = -18 00324 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 00325 INFO = -20 00326 END IF 00327 IF( INFO.NE.0 ) THEN 00328 CALL XERBLA( 'SGGSVP', -INFO ) 00329 RETURN 00330 END IF 00331 * 00332 * QR with column pivoting of B: B*P = V*( S11 S12 ) 00333 * ( 0 0 ) 00334 * 00335 DO 10 I = 1, N 00336 IWORK( I ) = 0 00337 10 CONTINUE 00338 CALL SGEQPF( P, N, B, LDB, IWORK, TAU, WORK, INFO ) 00339 * 00340 * Update A := A*P 00341 * 00342 CALL SLAPMT( FORWRD, M, N, A, LDA, IWORK ) 00343 * 00344 * Determine the effective rank of matrix B. 00345 * 00346 L = 0 00347 DO 20 I = 1, MIN( P, N ) 00348 IF( ABS( B( I, I ) ).GT.TOLB ) 00349 $ L = L + 1 00350 20 CONTINUE 00351 * 00352 IF( WANTV ) THEN 00353 * 00354 * Copy the details of V, and form V. 00355 * 00356 CALL SLASET( 'Full', P, P, ZERO, ZERO, V, LDV ) 00357 IF( P.GT.1 ) 00358 $ CALL SLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ), 00359 $ LDV ) 00360 CALL SORG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO ) 00361 END IF 00362 * 00363 * Clean up B 00364 * 00365 DO 40 J = 1, L - 1 00366 DO 30 I = J + 1, L 00367 B( I, J ) = ZERO 00368 30 CONTINUE 00369 40 CONTINUE 00370 IF( P.GT.L ) 00371 $ CALL SLASET( 'Full', P-L, N, ZERO, ZERO, B( L+1, 1 ), LDB ) 00372 * 00373 IF( WANTQ ) THEN 00374 * 00375 * Set Q = I and Update Q := Q*P 00376 * 00377 CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) 00378 CALL SLAPMT( FORWRD, N, N, Q, LDQ, IWORK ) 00379 END IF 00380 * 00381 IF( P.GE.L .AND. N.NE.L ) THEN 00382 * 00383 * RQ factorization of (S11 S12): ( S11 S12 ) = ( 0 S12 )*Z 00384 * 00385 CALL SGERQ2( L, N, B, LDB, TAU, WORK, INFO ) 00386 * 00387 * Update A := A*Z**T 00388 * 00389 CALL SORMR2( 'Right', 'Transpose', M, N, L, B, LDB, TAU, A, 00390 $ LDA, WORK, INFO ) 00391 * 00392 IF( WANTQ ) THEN 00393 * 00394 * Update Q := Q*Z**T 00395 * 00396 CALL SORMR2( 'Right', 'Transpose', N, N, L, B, LDB, TAU, Q, 00397 $ LDQ, WORK, INFO ) 00398 END IF 00399 * 00400 * Clean up B 00401 * 00402 CALL SLASET( 'Full', L, N-L, ZERO, ZERO, B, LDB ) 00403 DO 60 J = N - L + 1, N 00404 DO 50 I = J - N + L + 1, L 00405 B( I, J ) = ZERO 00406 50 CONTINUE 00407 60 CONTINUE 00408 * 00409 END IF 00410 * 00411 * Let N-L L 00412 * A = ( A11 A12 ) M, 00413 * 00414 * then the following does the complete QR decomposition of A11: 00415 * 00416 * A11 = U*( 0 T12 )*P1**T 00417 * ( 0 0 ) 00418 * 00419 DO 70 I = 1, N - L 00420 IWORK( I ) = 0 00421 70 CONTINUE 00422 CALL SGEQPF( M, N-L, A, LDA, IWORK, TAU, WORK, INFO ) 00423 * 00424 * Determine the effective rank of A11 00425 * 00426 K = 0 00427 DO 80 I = 1, MIN( M, N-L ) 00428 IF( ABS( A( I, I ) ).GT.TOLA ) 00429 $ K = K + 1 00430 80 CONTINUE 00431 * 00432 * Update A12 := U**T*A12, where A12 = A( 1:M, N-L+1:N ) 00433 * 00434 CALL SORM2R( 'Left', 'Transpose', M, L, MIN( M, N-L ), A, LDA, 00435 $ TAU, A( 1, N-L+1 ), LDA, WORK, INFO ) 00436 * 00437 IF( WANTU ) THEN 00438 * 00439 * Copy the details of U, and form U 00440 * 00441 CALL SLASET( 'Full', M, M, ZERO, ZERO, U, LDU ) 00442 IF( M.GT.1 ) 00443 $ CALL SLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ), 00444 $ LDU ) 00445 CALL SORG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO ) 00446 END IF 00447 * 00448 IF( WANTQ ) THEN 00449 * 00450 * Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1 00451 * 00452 CALL SLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK ) 00453 END IF 00454 * 00455 * Clean up A: set the strictly lower triangular part of 00456 * A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0. 00457 * 00458 DO 100 J = 1, K - 1 00459 DO 90 I = J + 1, K 00460 A( I, J ) = ZERO 00461 90 CONTINUE 00462 100 CONTINUE 00463 IF( M.GT.K ) 00464 $ CALL SLASET( 'Full', M-K, N-L, ZERO, ZERO, A( K+1, 1 ), LDA ) 00465 * 00466 IF( N-L.GT.K ) THEN 00467 * 00468 * RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1 00469 * 00470 CALL SGERQ2( K, N-L, A, LDA, TAU, WORK, INFO ) 00471 * 00472 IF( WANTQ ) THEN 00473 * 00474 * Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**T 00475 * 00476 CALL SORMR2( 'Right', 'Transpose', N, N-L, K, A, LDA, TAU, 00477 $ Q, LDQ, WORK, INFO ) 00478 END IF 00479 * 00480 * Clean up A 00481 * 00482 CALL SLASET( 'Full', K, N-L-K, ZERO, ZERO, A, LDA ) 00483 DO 120 J = N - L - K + 1, N - L 00484 DO 110 I = J - N + L + K + 1, K 00485 A( I, J ) = ZERO 00486 110 CONTINUE 00487 120 CONTINUE 00488 * 00489 END IF 00490 * 00491 IF( M.GT.K ) THEN 00492 * 00493 * QR factorization of A( K+1:M,N-L+1:N ) 00494 * 00495 CALL SGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO ) 00496 * 00497 IF( WANTU ) THEN 00498 * 00499 * Update U(:,K+1:M) := U(:,K+1:M)*U1 00500 * 00501 CALL SORM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ), 00502 $ A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU, 00503 $ WORK, INFO ) 00504 END IF 00505 * 00506 * Clean up 00507 * 00508 DO 140 J = N - L + 1, N 00509 DO 130 I = J - N + K + L + 1, M 00510 A( I, J ) = ZERO 00511 130 CONTINUE 00512 140 CONTINUE 00513 * 00514 END IF 00515 * 00516 RETURN 00517 * 00518 * End of SGGSVP 00519 * 00520 END