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