 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, 00002 $ TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, 00003 $ IWORK, TAU, WORK, INFO ) 00004 * 00005 * -- LAPACK routine (version 3.3.1) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * -- April 2011 -- 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER JOBQ, JOBU, JOBV 00012 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P 00013 REAL TOLA, TOLB 00014 * .. 00015 * .. Array Arguments .. 00016 INTEGER IWORK( * ) 00017 REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00018 $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * SGGSVP computes orthogonal matrices U, V and Q such that 00025 * 00026 * N-K-L K L 00027 * U**T*A*Q = K ( 0 A12 A13 ) if M-K-L >= 0; 00028 * L ( 0 0 A23 ) 00029 * M-K-L ( 0 0 0 ) 00030 * 00031 * N-K-L K L 00032 * = K ( 0 A12 A13 ) if M-K-L < 0; 00033 * M-K ( 0 0 A23 ) 00034 * 00035 * N-K-L K L 00036 * V**T*B*Q = L ( 0 0 B13 ) 00037 * P-L ( 0 0 0 ) 00038 * 00039 * where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular 00040 * upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, 00041 * otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective 00042 * numerical rank of the (M+P)-by-N matrix (A**T,B**T)**T. 00043 * 00044 * This decomposition is the preprocessing step for computing the 00045 * Generalized Singular Value Decomposition (GSVD), see subroutine 00046 * SGGSVD. 00047 * 00048 * Arguments 00049 * ========= 00050 * 00051 * JOBU (input) CHARACTER*1 00052 * = 'U': Orthogonal matrix U is computed; 00053 * = 'N': U is not computed. 00054 * 00055 * JOBV (input) CHARACTER*1 00056 * = 'V': Orthogonal matrix V is computed; 00057 * = 'N': V is not computed. 00058 * 00059 * JOBQ (input) CHARACTER*1 00060 * = 'Q': Orthogonal matrix Q is computed; 00061 * = 'N': Q is not computed. 00062 * 00063 * M (input) INTEGER 00064 * The number of rows of the matrix A. M >= 0. 00065 * 00066 * P (input) INTEGER 00067 * The number of rows of the matrix B. P >= 0. 00068 * 00069 * N (input) INTEGER 00070 * The number of columns of the matrices A and B. N >= 0. 00071 * 00072 * A (input/output) REAL array, dimension (LDA,N) 00073 * On entry, the M-by-N matrix A. 00074 * On exit, A contains the triangular (or trapezoidal) matrix 00075 * described in the Purpose section. 00076 * 00077 * LDA (input) INTEGER 00078 * The leading dimension of the array A. LDA >= max(1,M). 00079 * 00080 * B (input/output) REAL array, dimension (LDB,N) 00081 * On entry, the P-by-N matrix B. 00082 * On exit, B contains the triangular matrix described in 00083 * the Purpose section. 00084 * 00085 * LDB (input) INTEGER 00086 * The leading dimension of the array B. LDB >= max(1,P). 00087 * 00088 * TOLA (input) REAL 00089 * TOLB (input) REAL 00090 * TOLA and TOLB are the thresholds to determine the effective 00091 * numerical rank of matrix B and a subblock of A. Generally, 00092 * they are set to 00093 * TOLA = MAX(M,N)*norm(A)*MACHEPS, 00094 * TOLB = MAX(P,N)*norm(B)*MACHEPS. 00095 * The size of TOLA and TOLB may affect the size of backward 00096 * errors of the decomposition. 00097 * 00098 * K (output) INTEGER 00099 * L (output) INTEGER 00100 * On exit, K and L specify the dimension of the subblocks 00101 * described in Purpose section. 00102 * K + L = effective numerical rank of (A**T,B**T)**T. 00103 * 00104 * U (output) REAL array, dimension (LDU,M) 00105 * If JOBU = 'U', U contains the orthogonal matrix U. 00106 * If JOBU = 'N', U is not referenced. 00107 * 00108 * LDU (input) INTEGER 00109 * The leading dimension of the array U. LDU >= max(1,M) if 00110 * JOBU = 'U'; LDU >= 1 otherwise. 00111 * 00112 * V (output) REAL array, dimension (LDV,P) 00113 * If JOBV = 'V', V contains the orthogonal matrix V. 00114 * If JOBV = 'N', V is not referenced. 00115 * 00116 * LDV (input) INTEGER 00117 * The leading dimension of the array V. LDV >= max(1,P) if 00118 * JOBV = 'V'; LDV >= 1 otherwise. 00119 * 00120 * Q (output) REAL array, dimension (LDQ,N) 00121 * If JOBQ = 'Q', Q contains the orthogonal matrix Q. 00122 * If JOBQ = 'N', Q is not referenced. 00123 * 00124 * LDQ (input) INTEGER 00125 * The leading dimension of the array Q. LDQ >= max(1,N) if 00126 * JOBQ = 'Q'; LDQ >= 1 otherwise. 00127 * 00128 * IWORK (workspace) INTEGER array, dimension (N) 00129 * 00130 * TAU (workspace) REAL array, dimension (N) 00131 * 00132 * WORK (workspace) REAL array, dimension (max(3*N,M,P)) 00133 * 00134 * INFO (output) INTEGER 00135 * = 0: successful exit 00136 * < 0: if INFO = -i, the i-th argument had an illegal value. 00137 * 00138 * 00139 * Further Details 00140 * =============== 00141 * 00142 * The subroutine uses LAPACK subroutine SGEQPF for the QR factorization 00143 * with column pivoting to detect the effective numerical rank of the 00144 * a matrix. It may be replaced by a better rank determination strategy. 00145 * 00146 * ===================================================================== 00147 * 00148 * .. Parameters .. 00149 REAL ZERO, ONE 00150 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00151 * .. 00152 * .. Local Scalars .. 00153 LOGICAL FORWRD, WANTQ, WANTU, WANTV 00154 INTEGER I, J 00155 * .. 00156 * .. External Functions .. 00157 LOGICAL LSAME 00158 EXTERNAL LSAME 00159 * .. 00160 * .. External Subroutines .. 00161 EXTERNAL SGEQPF, SGEQR2, SGERQ2, SLACPY, SLAPMT, SLASET, 00162 $ SORG2R, SORM2R, SORMR2, XERBLA 00163 * .. 00164 * .. Intrinsic Functions .. 00165 INTRINSIC ABS, MAX, MIN 00166 * .. 00167 * .. Executable Statements .. 00168 * 00169 * Test the input parameters 00170 * 00171 WANTU = LSAME( JOBU, 'U' ) 00172 WANTV = LSAME( JOBV, 'V' ) 00173 WANTQ = LSAME( JOBQ, 'Q' ) 00174 FORWRD = .TRUE. 00175 * 00176 INFO = 0 00177 IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN 00178 INFO = -1 00179 ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN 00180 INFO = -2 00181 ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN 00182 INFO = -3 00183 ELSE IF( M.LT.0 ) THEN 00184 INFO = -4 00185 ELSE IF( P.LT.0 ) THEN 00186 INFO = -5 00187 ELSE IF( N.LT.0 ) THEN 00188 INFO = -6 00189 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00190 INFO = -8 00191 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 00192 INFO = -10 00193 ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN 00194 INFO = -16 00195 ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN 00196 INFO = -18 00197 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 00198 INFO = -20 00199 END IF 00200 IF( INFO.NE.0 ) THEN 00201 CALL XERBLA( 'SGGSVP', -INFO ) 00202 RETURN 00203 END IF 00204 * 00205 * QR with column pivoting of B: B*P = V*( S11 S12 ) 00206 * ( 0 0 ) 00207 * 00208 DO 10 I = 1, N 00209 IWORK( I ) = 0 00210 10 CONTINUE 00211 CALL SGEQPF( P, N, B, LDB, IWORK, TAU, WORK, INFO ) 00212 * 00213 * Update A := A*P 00214 * 00215 CALL SLAPMT( FORWRD, M, N, A, LDA, IWORK ) 00216 * 00217 * Determine the effective rank of matrix B. 00218 * 00219 L = 0 00220 DO 20 I = 1, MIN( P, N ) 00221 IF( ABS( B( I, I ) ).GT.TOLB ) 00222 $ L = L + 1 00223 20 CONTINUE 00224 * 00225 IF( WANTV ) THEN 00226 * 00227 * Copy the details of V, and form V. 00228 * 00229 CALL SLASET( 'Full', P, P, ZERO, ZERO, V, LDV ) 00230 IF( P.GT.1 ) 00231 $ CALL SLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ), 00232 $ LDV ) 00233 CALL SORG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO ) 00234 END IF 00235 * 00236 * Clean up B 00237 * 00238 DO 40 J = 1, L - 1 00239 DO 30 I = J + 1, L 00240 B( I, J ) = ZERO 00241 30 CONTINUE 00242 40 CONTINUE 00243 IF( P.GT.L ) 00244 $ CALL SLASET( 'Full', P-L, N, ZERO, ZERO, B( L+1, 1 ), LDB ) 00245 * 00246 IF( WANTQ ) THEN 00247 * 00248 * Set Q = I and Update Q := Q*P 00249 * 00250 CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) 00251 CALL SLAPMT( FORWRD, N, N, Q, LDQ, IWORK ) 00252 END IF 00253 * 00254 IF( P.GE.L .AND. N.NE.L ) THEN 00255 * 00256 * RQ factorization of (S11 S12): ( S11 S12 ) = ( 0 S12 )*Z 00257 * 00258 CALL SGERQ2( L, N, B, LDB, TAU, WORK, INFO ) 00259 * 00260 * Update A := A*Z**T 00261 * 00262 CALL SORMR2( 'Right', 'Transpose', M, N, L, B, LDB, TAU, A, 00263 $ LDA, WORK, INFO ) 00264 * 00265 IF( WANTQ ) THEN 00266 * 00267 * Update Q := Q*Z**T 00268 * 00269 CALL SORMR2( 'Right', 'Transpose', N, N, L, B, LDB, TAU, Q, 00270 $ LDQ, WORK, INFO ) 00271 END IF 00272 * 00273 * Clean up B 00274 * 00275 CALL SLASET( 'Full', L, N-L, ZERO, ZERO, B, LDB ) 00276 DO 60 J = N - L + 1, N 00277 DO 50 I = J - N + L + 1, L 00278 B( I, J ) = ZERO 00279 50 CONTINUE 00280 60 CONTINUE 00281 * 00282 END IF 00283 * 00284 * Let N-L L 00285 * A = ( A11 A12 ) M, 00286 * 00287 * then the following does the complete QR decomposition of A11: 00288 * 00289 * A11 = U*( 0 T12 )*P1**T 00290 * ( 0 0 ) 00291 * 00292 DO 70 I = 1, N - L 00293 IWORK( I ) = 0 00294 70 CONTINUE 00295 CALL SGEQPF( M, N-L, A, LDA, IWORK, TAU, WORK, INFO ) 00296 * 00297 * Determine the effective rank of A11 00298 * 00299 K = 0 00300 DO 80 I = 1, MIN( M, N-L ) 00301 IF( ABS( A( I, I ) ).GT.TOLA ) 00302 $ K = K + 1 00303 80 CONTINUE 00304 * 00305 * Update A12 := U**T*A12, where A12 = A( 1:M, N-L+1:N ) 00306 * 00307 CALL SORM2R( 'Left', 'Transpose', M, L, MIN( M, N-L ), A, LDA, 00308 $ TAU, A( 1, N-L+1 ), LDA, WORK, INFO ) 00309 * 00310 IF( WANTU ) THEN 00311 * 00312 * Copy the details of U, and form U 00313 * 00314 CALL SLASET( 'Full', M, M, ZERO, ZERO, U, LDU ) 00315 IF( M.GT.1 ) 00316 $ CALL SLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ), 00317 $ LDU ) 00318 CALL SORG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO ) 00319 END IF 00320 * 00321 IF( WANTQ ) THEN 00322 * 00323 * Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1 00324 * 00325 CALL SLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK ) 00326 END IF 00327 * 00328 * Clean up A: set the strictly lower triangular part of 00329 * A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0. 00330 * 00331 DO 100 J = 1, K - 1 00332 DO 90 I = J + 1, K 00333 A( I, J ) = ZERO 00334 90 CONTINUE 00335 100 CONTINUE 00336 IF( M.GT.K ) 00337 $ CALL SLASET( 'Full', M-K, N-L, ZERO, ZERO, A( K+1, 1 ), LDA ) 00338 * 00339 IF( N-L.GT.K ) THEN 00340 * 00341 * RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1 00342 * 00343 CALL SGERQ2( K, N-L, A, LDA, TAU, WORK, INFO ) 00344 * 00345 IF( WANTQ ) THEN 00346 * 00347 * Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**T 00348 * 00349 CALL SORMR2( 'Right', 'Transpose', N, N-L, K, A, LDA, TAU, 00350 $ Q, LDQ, WORK, INFO ) 00351 END IF 00352 * 00353 * Clean up A 00354 * 00355 CALL SLASET( 'Full', K, N-L-K, ZERO, ZERO, A, LDA ) 00356 DO 120 J = N - L - K + 1, N - L 00357 DO 110 I = J - N + L + K + 1, K 00358 A( I, J ) = ZERO 00359 110 CONTINUE 00360 120 CONTINUE 00361 * 00362 END IF 00363 * 00364 IF( M.GT.K ) THEN 00365 * 00366 * QR factorization of A( K+1:M,N-L+1:N ) 00367 * 00368 CALL SGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO ) 00369 * 00370 IF( WANTU ) THEN 00371 * 00372 * Update U(:,K+1:M) := U(:,K+1:M)*U1 00373 * 00374 CALL SORM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ), 00375 $ A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU, 00376 $ WORK, INFO ) 00377 END IF 00378 * 00379 * Clean up 00380 * 00381 DO 140 J = N - L + 1, N 00382 DO 130 I = J - N + K + L + 1, M 00383 A( I, J ) = ZERO 00384 130 CONTINUE 00385 140 CONTINUE 00386 * 00387 END IF 00388 * 00389 RETURN 00390 * 00391 * End of SGGSVP 00392 * 00393 END
1.7.3 
