3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download SGGSVP + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sggsvp.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sggsvp.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sggsvp.f">
21 * SUBROUTINE SGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
22 * TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
23 * IWORK, TAU, WORK, INFO )
25 * .. Scalar Arguments ..
26 * CHARACTER JOBQ, JOBU, JOBV
27 * INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
30 * .. Array Arguments ..
32 * REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
33 * $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
42 *> This routine is deprecated and has been replaced by routine SGGSVP3.
44 *> SGGSVP computes orthogonal matrices U, V and Q such that
47 *> U**T*A*Q = K ( 0 A12 A13 ) if M-K-L >= 0;
52 *> = K ( 0 A12 A13 ) if M-K-L < 0;
56 *> V**T*B*Q = L ( 0 0 B13 )
59 *> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
60 *> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
61 *> otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective
62 *> numerical rank of the (M+P)-by-N matrix (A**T,B**T)**T.
64 *> This decomposition is the preprocessing step for computing the
65 *> Generalized Singular Value Decomposition (GSVD), see subroutine
74 *> JOBU is CHARACTER*1
75 *> = 'U': Orthogonal matrix U is computed;
76 *> = 'N': U is not computed.
81 *> JOBV is CHARACTER*1
82 *> = 'V': Orthogonal matrix V is computed;
83 *> = 'N': V is not computed.
88 *> JOBQ is CHARACTER*1
89 *> = 'Q': Orthogonal matrix Q is computed;
90 *> = 'N': Q is not computed.
96 *> The number of rows of the matrix A. M >= 0.
102 *> The number of rows of the matrix B. P >= 0.
108 *> The number of columns of the matrices A and B. N >= 0.
113 *> A is REAL array, dimension (LDA,N)
114 *> On entry, the M-by-N matrix A.
115 *> On exit, A contains the triangular (or trapezoidal) matrix
116 *> described in the Purpose section.
122 *> The leading dimension of the array A. LDA >= max(1,M).
127 *> B is REAL array, dimension (LDB,N)
128 *> On entry, the P-by-N matrix B.
129 *> On exit, B contains the triangular matrix described in
130 *> the Purpose section.
136 *> The leading dimension of the array B. LDB >= max(1,P).
148 *> TOLA and TOLB are the thresholds to determine the effective
149 *> numerical rank of matrix B and a subblock of A. Generally,
151 *> TOLA = MAX(M,N)*norm(A)*MACHEPS,
152 *> TOLB = MAX(P,N)*norm(B)*MACHEPS.
153 *> The size of TOLA and TOLB may affect the size of backward
154 *> errors of the decomposition.
166 *> On exit, K and L specify the dimension of the subblocks
167 *> described in Purpose section.
168 *> K + L = effective numerical rank of (A**T,B**T)**T.
173 *> U is REAL array, dimension (LDU,M)
174 *> If JOBU = 'U', U contains the orthogonal matrix U.
175 *> If JOBU = 'N', U is not referenced.
181 *> The leading dimension of the array U. LDU >= max(1,M) if
182 *> JOBU = 'U'; LDU >= 1 otherwise.
187 *> V is REAL array, dimension (LDV,P)
188 *> If JOBV = 'V', V contains the orthogonal matrix V.
189 *> If JOBV = 'N', V is not referenced.
195 *> The leading dimension of the array V. LDV >= max(1,P) if
196 *> JOBV = 'V'; LDV >= 1 otherwise.
201 *> Q is REAL array, dimension (LDQ,N)
202 *> If JOBQ = 'Q', Q contains the orthogonal matrix Q.
203 *> If JOBQ = 'N', Q is not referenced.
209 *> The leading dimension of the array Q. LDQ >= max(1,N) if
210 *> JOBQ = 'Q'; LDQ >= 1 otherwise.
215 *> IWORK is INTEGER array, dimension (N)
220 *> TAU is REAL array, dimension (N)
225 *> WORK is REAL array, dimension (max(3*N,M,P))
231 *> = 0: successful exit
232 *> < 0: if INFO = -i, the i-th argument had an illegal value.
238 *> \author Univ. of Tennessee
239 *> \author Univ. of California Berkeley
240 *> \author Univ. of Colorado Denver
243 *> \date November 2011
245 *> \ingroup realOTHERcomputational
247 *> \par Further Details:
248 * =====================
250 *> The subroutine uses LAPACK subroutine SGEQPF for the QR factorization
251 *> with column pivoting to detect the effective numerical rank of the
252 *> a matrix. It may be replaced by a better rank determination strategy.
254 * =====================================================================
255 SUBROUTINE SGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
256 $ TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
257 $ IWORK, TAU, WORK, INFO )
259 * -- LAPACK computational routine (version 3.4.0) --
260 * -- LAPACK is a software package provided by Univ. of Tennessee, --
261 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
264 * .. Scalar Arguments ..
265 CHARACTER JOBQ, JOBU, JOBV
266 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
269 * .. Array Arguments ..
271 REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
272 $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
275 * =====================================================================
279 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
281 * .. Local Scalars ..
282 LOGICAL FORWRD, WANTQ, WANTU, WANTV
285 * .. External Functions ..
289 * .. External Subroutines ..
290 EXTERNAL SGEQPF, SGEQR2, SGERQ2, SLACPY, SLAPMT, SLASET,
291 $ SORG2R, SORM2R, SORMR2, XERBLA
293 * .. Intrinsic Functions ..
294 INTRINSIC ABS, MAX, MIN
296 * .. Executable Statements ..
298 * Test the input parameters
300 WANTU = LSAME( JOBU, 'U' )
301 WANTV = LSAME( JOBV, 'V' )
302 WANTQ = LSAME( JOBQ, 'Q' )
306 IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN
308 ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN
310 ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN
312 ELSE IF( M.LT.0 ) THEN
314 ELSE IF( P.LT.0 ) THEN
316 ELSE IF( N.LT.0 ) THEN
318 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
320 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
322 ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
324 ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
326 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
330 CALL XERBLA( 'SGGSVP', -INFO )
334 * QR with column pivoting of B: B*P = V*( S11 S12 )
340 CALL SGEQPF( P, N, B, LDB, IWORK, TAU, WORK, INFO )
344 CALL SLAPMT( FORWRD, M, N, A, LDA, IWORK )
346 * Determine the effective rank of matrix B.
349 DO 20 I = 1, MIN( P, N )
350 IF( ABS( B( I, I ) ).GT.TOLB )
356 * Copy the details of V, and form V.
358 CALL SLASET( 'Full', P, P, ZERO, ZERO, V, LDV )
360 $ CALL SLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ),
362 CALL SORG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO )
373 $ CALL SLASET( 'Full', P-L, N, ZERO, ZERO, B( L+1, 1 ), LDB )
377 * Set Q = I and Update Q := Q*P
379 CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
380 CALL SLAPMT( FORWRD, N, N, Q, LDQ, IWORK )
383 IF( P.GE.L .AND. N.NE.L ) THEN
385 * RQ factorization of (S11 S12): ( S11 S12 ) = ( 0 S12 )*Z
387 CALL SGERQ2( L, N, B, LDB, TAU, WORK, INFO )
391 CALL SORMR2( 'Right', 'Transpose', M, N, L, B, LDB, TAU, A,
398 CALL SORMR2( 'Right', 'Transpose', N, N, L, B, LDB, TAU, Q,
404 CALL SLASET( 'Full', L, N-L, ZERO, ZERO, B, LDB )
405 DO 60 J = N - L + 1, N
406 DO 50 I = J - N + L + 1, L
416 * then the following does the complete QR decomposition of A11:
418 * A11 = U*( 0 T12 )*P1**T
424 CALL SGEQPF( M, N-L, A, LDA, IWORK, TAU, WORK, INFO )
426 * Determine the effective rank of A11
429 DO 80 I = 1, MIN( M, N-L )
430 IF( ABS( A( I, I ) ).GT.TOLA )
434 * Update A12 := U**T*A12, where A12 = A( 1:M, N-L+1:N )
436 CALL SORM2R( 'Left', 'Transpose', M, L, MIN( M, N-L ), A, LDA,
437 $ TAU, A( 1, N-L+1 ), LDA, WORK, INFO )
441 * Copy the details of U, and form U
443 CALL SLASET( 'Full', M, M, ZERO, ZERO, U, LDU )
445 $ CALL SLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ),
447 CALL SORG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO )
452 * Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1
454 CALL SLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK )
457 * Clean up A: set the strictly lower triangular part of
458 * A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0.
466 $ CALL SLASET( 'Full', M-K, N-L, ZERO, ZERO, A( K+1, 1 ), LDA )
470 * RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1
472 CALL SGERQ2( K, N-L, A, LDA, TAU, WORK, INFO )
476 * Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**T
478 CALL SORMR2( 'Right', 'Transpose', N, N-L, K, A, LDA, TAU,
479 $ Q, LDQ, WORK, INFO )
484 CALL SLASET( 'Full', K, N-L-K, ZERO, ZERO, A, LDA )
485 DO 120 J = N - L - K + 1, N - L
486 DO 110 I = J - N + L + K + 1, K
495 * QR factorization of A( K+1:M,N-L+1:N )
497 CALL SGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO )
501 * Update U(:,K+1:M) := U(:,K+1:M)*U1
503 CALL SORM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ),
504 $ A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU,
510 DO 140 J = N - L + 1, N
511 DO 130 I = J - N + K + L + 1, M