LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cggsvp3.f
Go to the documentation of this file.
1*> \brief \b CGGSVP3
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CGGSVP3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggsvp3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggsvp3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggsvp3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
22* TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
23* IWORK, RWORK, TAU, WORK, LWORK, INFO )
24*
25* .. Scalar Arguments ..
26* CHARACTER JOBQ, JOBU, JOBV
27* INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, LWORK
28* REAL TOLA, TOLB
29* ..
30* .. Array Arguments ..
31* INTEGER IWORK( * )
32* REAL RWORK( * )
33* COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
34* $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
35* ..
36*
37*
38*> \par Purpose:
39* =============
40*>
41*> \verbatim
42*>
43*> CGGSVP3 computes unitary matrices U, V and Q such that
44*>
45*> N-K-L K L
46*> U**H*A*Q = K ( 0 A12 A13 ) if M-K-L >= 0;
47*> L ( 0 0 A23 )
48*> M-K-L ( 0 0 0 )
49*>
50*> N-K-L K L
51*> = K ( 0 A12 A13 ) if M-K-L < 0;
52*> M-K ( 0 0 A23 )
53*>
54*> N-K-L K L
55*> V**H*B*Q = L ( 0 0 B13 )
56*> P-L ( 0 0 0 )
57*>
58*> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
59*> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
60*> otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective
61*> numerical rank of the (M+P)-by-N matrix (A**H,B**H)**H.
62*>
63*> This decomposition is the preprocessing step for computing the
64*> Generalized Singular Value Decomposition (GSVD), see subroutine
65*> CGGSVD3.
66*> \endverbatim
67*
68* Arguments:
69* ==========
70*
71*> \param[in] JOBU
72*> \verbatim
73*> JOBU is CHARACTER*1
74*> = 'U': Unitary matrix U is computed;
75*> = 'N': U is not computed.
76*> \endverbatim
77*>
78*> \param[in] JOBV
79*> \verbatim
80*> JOBV is CHARACTER*1
81*> = 'V': Unitary matrix V is computed;
82*> = 'N': V is not computed.
83*> \endverbatim
84*>
85*> \param[in] JOBQ
86*> \verbatim
87*> JOBQ is CHARACTER*1
88*> = 'Q': Unitary matrix Q is computed;
89*> = 'N': Q is not computed.
90*> \endverbatim
91*>
92*> \param[in] M
93*> \verbatim
94*> M is INTEGER
95*> The number of rows of the matrix A. M >= 0.
96*> \endverbatim
97*>
98*> \param[in] P
99*> \verbatim
100*> P is INTEGER
101*> The number of rows of the matrix B. P >= 0.
102*> \endverbatim
103*>
104*> \param[in] N
105*> \verbatim
106*> N is INTEGER
107*> The number of columns of the matrices A and B. N >= 0.
108*> \endverbatim
109*>
110*> \param[in,out] A
111*> \verbatim
112*> A is COMPLEX array, dimension (LDA,N)
113*> On entry, the M-by-N matrix A.
114*> On exit, A contains the triangular (or trapezoidal) matrix
115*> described in the Purpose section.
116*> \endverbatim
117*>
118*> \param[in] LDA
119*> \verbatim
120*> LDA is INTEGER
121*> The leading dimension of the array A. LDA >= max(1,M).
122*> \endverbatim
123*>
124*> \param[in,out] B
125*> \verbatim
126*> B is COMPLEX array, dimension (LDB,N)
127*> On entry, the P-by-N matrix B.
128*> On exit, B contains the triangular matrix described in
129*> the Purpose section.
130*> \endverbatim
131*>
132*> \param[in] LDB
133*> \verbatim
134*> LDB is INTEGER
135*> The leading dimension of the array B. LDB >= max(1,P).
136*> \endverbatim
137*>
138*> \param[in] TOLA
139*> \verbatim
140*> TOLA is REAL
141*> \endverbatim
142*>
143*> \param[in] TOLB
144*> \verbatim
145*> TOLB is REAL
146*>
147*> TOLA and TOLB are the thresholds to determine the effective
148*> numerical rank of matrix B and a subblock of A. Generally,
149*> they are set to
150*> TOLA = MAX(M,N)*norm(A)*MACHEPS,
151*> TOLB = MAX(P,N)*norm(B)*MACHEPS.
152*> The size of TOLA and TOLB may affect the size of backward
153*> errors of the decomposition.
154*> \endverbatim
155*>
156*> \param[out] K
157*> \verbatim
158*> K is INTEGER
159*> \endverbatim
160*>
161*> \param[out] L
162*> \verbatim
163*> L is INTEGER
164*>
165*> On exit, K and L specify the dimension of the subblocks
166*> described in Purpose section.
167*> K + L = effective numerical rank of (A**H,B**H)**H.
168*> \endverbatim
169*>
170*> \param[out] U
171*> \verbatim
172*> U is COMPLEX array, dimension (LDU,M)
173*> If JOBU = 'U', U contains the unitary matrix U.
174*> If JOBU = 'N', U is not referenced.
175*> \endverbatim
176*>
177*> \param[in] LDU
178*> \verbatim
179*> LDU is INTEGER
180*> The leading dimension of the array U. LDU >= max(1,M) if
181*> JOBU = 'U'; LDU >= 1 otherwise.
182*> \endverbatim
183*>
184*> \param[out] V
185*> \verbatim
186*> V is COMPLEX array, dimension (LDV,P)
187*> If JOBV = 'V', V contains the unitary matrix V.
188*> If JOBV = 'N', V is not referenced.
189*> \endverbatim
190*>
191*> \param[in] LDV
192*> \verbatim
193*> LDV is INTEGER
194*> The leading dimension of the array V. LDV >= max(1,P) if
195*> JOBV = 'V'; LDV >= 1 otherwise.
196*> \endverbatim
197*>
198*> \param[out] Q
199*> \verbatim
200*> Q is COMPLEX array, dimension (LDQ,N)
201*> If JOBQ = 'Q', Q contains the unitary matrix Q.
202*> If JOBQ = 'N', Q is not referenced.
203*> \endverbatim
204*>
205*> \param[in] LDQ
206*> \verbatim
207*> LDQ is INTEGER
208*> The leading dimension of the array Q. LDQ >= max(1,N) if
209*> JOBQ = 'Q'; LDQ >= 1 otherwise.
210*> \endverbatim
211*>
212*> \param[out] IWORK
213*> \verbatim
214*> IWORK is INTEGER array, dimension (N)
215*> \endverbatim
216*>
217*> \param[out] RWORK
218*> \verbatim
219*> RWORK is REAL array, dimension (2*N)
220*> \endverbatim
221*>
222*> \param[out] TAU
223*> \verbatim
224*> TAU is COMPLEX array, dimension (N)
225*> \endverbatim
226*>
227*> \param[out] WORK
228*> \verbatim
229*> WORK is COMPLEX array, dimension (MAX(1,LWORK))
230*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
231*> \endverbatim
232*>
233*> \param[in] LWORK
234*> \verbatim
235*> LWORK is INTEGER
236*> The dimension of the array WORK.
237*>
238*> If LWORK = -1, then a workspace query is assumed; the routine
239*> only calculates the optimal size of the WORK array, returns
240*> this value as the first entry of the WORK array, and no error
241*> message related to LWORK is issued by XERBLA.
242*> \endverbatim
243*>
244*> \param[out] INFO
245*> \verbatim
246*> INFO is INTEGER
247*> = 0: successful exit
248*> < 0: if INFO = -i, the i-th argument had an illegal value.
249*> \endverbatim
250*
251* Authors:
252* ========
253*
254*> \author Univ. of Tennessee
255*> \author Univ. of California Berkeley
256*> \author Univ. of Colorado Denver
257*> \author NAG Ltd.
258*
259*> \ingroup ggsvp3
260*
261*> \par Further Details:
262* =====================
263*>
264*> \verbatim
265*>
266*> The subroutine uses LAPACK subroutine CGEQP3 for the QR factorization
267*> with column pivoting to detect the effective numerical rank of the
268*> a matrix. It may be replaced by a better rank determination strategy.
269*>
270*> CGGSVP3 replaces the deprecated subroutine CGGSVP.
271*>
272*> \endverbatim
273*>
274* =====================================================================
275 SUBROUTINE cggsvp3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
276 $ TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
277 $ IWORK, RWORK, TAU, WORK, LWORK, INFO )
278*
279* -- LAPACK computational routine --
280* -- LAPACK is a software package provided by Univ. of Tennessee, --
281* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
282*
283 IMPLICIT NONE
284*
285* .. Scalar Arguments ..
286 CHARACTER JOBQ, JOBU, JOBV
287 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P,
288 $ lwork
289 REAL TOLA, TOLB
290* ..
291* .. Array Arguments ..
292 INTEGER IWORK( * )
293 REAL RWORK( * )
294 COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
295 $ tau( * ), u( ldu, * ), v( ldv, * ), work( * )
296* ..
297*
298* =====================================================================
299*
300* .. Parameters ..
301 COMPLEX CZERO, CONE
302 PARAMETER ( CZERO = ( 0.0e+0, 0.0e+0 ),
303 $ cone = ( 1.0e+0, 0.0e+0 ) )
304* ..
305* .. Local Scalars ..
306 LOGICAL FORWRD, WANTQ, WANTU, WANTV, LQUERY
307 INTEGER I, J, LWKOPT
308* ..
309* .. External Functions ..
310 LOGICAL LSAME
311 EXTERNAL LSAME
312* ..
313* .. External Subroutines ..
314 EXTERNAL cgeqp3, cgeqr2, cgerq2, clacpy, clapmt,
316* ..
317* .. Intrinsic Functions ..
318 INTRINSIC abs, aimag, max, min, real
319* ..
320* .. Executable Statements ..
321*
322* Test the input parameters
323*
324 wantu = lsame( jobu, 'U' )
325 wantv = lsame( jobv, 'V' )
326 wantq = lsame( jobq, 'Q' )
327 forwrd = .true.
328 lquery = ( lwork.EQ.-1 )
329 lwkopt = 1
330*
331* Test the input arguments
332*
333 info = 0
334 IF( .NOT.( wantu .OR. lsame( jobu, 'N' ) ) ) THEN
335 info = -1
336 ELSE IF( .NOT.( wantv .OR. lsame( jobv, 'N' ) ) ) THEN
337 info = -2
338 ELSE IF( .NOT.( wantq .OR. lsame( jobq, 'N' ) ) ) THEN
339 info = -3
340 ELSE IF( m.LT.0 ) THEN
341 info = -4
342 ELSE IF( p.LT.0 ) THEN
343 info = -5
344 ELSE IF( n.LT.0 ) THEN
345 info = -6
346 ELSE IF( lda.LT.max( 1, m ) ) THEN
347 info = -8
348 ELSE IF( ldb.LT.max( 1, p ) ) THEN
349 info = -10
350 ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) ) THEN
351 info = -16
352 ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) ) THEN
353 info = -18
354 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
355 info = -20
356 ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
357 info = -24
358 END IF
359*
360* Compute workspace
361*
362 IF( info.EQ.0 ) THEN
363 CALL cgeqp3( p, n, b, ldb, iwork, tau, work, -1, rwork, info )
364 lwkopt = int( work( 1 ) )
365 IF( wantv ) THEN
366 lwkopt = max( lwkopt, p )
367 END IF
368 lwkopt = max( lwkopt, min( n, p ) )
369 lwkopt = max( lwkopt, m )
370 IF( wantq ) THEN
371 lwkopt = max( lwkopt, n )
372 END IF
373 CALL cgeqp3( m, n, a, lda, iwork, tau, work, -1, rwork, info )
374 lwkopt = max( lwkopt, int( work( 1 ) ) )
375 lwkopt = max( 1, lwkopt )
376 work( 1 ) = cmplx( lwkopt )
377 END IF
378*
379 IF( info.NE.0 ) THEN
380 CALL xerbla( 'CGGSVP3', -info )
381 RETURN
382 END IF
383 IF( lquery ) THEN
384 RETURN
385 ENDIF
386*
387* QR with column pivoting of B: B*P = V*( S11 S12 )
388* ( 0 0 )
389*
390 DO 10 i = 1, n
391 iwork( i ) = 0
392 10 CONTINUE
393 CALL cgeqp3( p, n, b, ldb, iwork, tau, work, lwork, rwork, info )
394*
395* Update A := A*P
396*
397 CALL clapmt( forwrd, m, n, a, lda, iwork )
398*
399* Determine the effective rank of matrix B.
400*
401 l = 0
402 DO 20 i = 1, min( p, n )
403 IF( abs( b( i, i ) ).GT.tolb )
404 $ l = l + 1
405 20 CONTINUE
406*
407 IF( wantv ) THEN
408*
409* Copy the details of V, and form V.
410*
411 CALL claset( 'Full', p, p, czero, czero, v, ldv )
412 IF( p.GT.1 )
413 $ CALL clacpy( 'Lower', p-1, n, b( 2, 1 ), ldb, v( 2, 1 ),
414 $ ldv )
415 CALL cung2r( p, p, min( p, n ), v, ldv, tau, work, info )
416 END IF
417*
418* Clean up B
419*
420 DO 40 j = 1, l - 1
421 DO 30 i = j + 1, l
422 b( i, j ) = czero
423 30 CONTINUE
424 40 CONTINUE
425 IF( p.GT.l )
426 $ CALL claset( 'Full', p-l, n, czero, czero, b( l+1, 1 ), ldb )
427*
428 IF( wantq ) THEN
429*
430* Set Q = I and Update Q := Q*P
431*
432 CALL claset( 'Full', n, n, czero, cone, q, ldq )
433 CALL clapmt( forwrd, n, n, q, ldq, iwork )
434 END IF
435*
436 IF( p.GE.l .AND. n.NE.l ) THEN
437*
438* RQ factorization of ( S11 S12 ) = ( 0 S12 )*Z
439*
440 CALL cgerq2( l, n, b, ldb, tau, work, info )
441*
442* Update A := A*Z**H
443*
444 CALL cunmr2( 'Right', 'Conjugate transpose', m, n, l, b, ldb,
445 $ tau, a, lda, work, info )
446 IF( wantq ) THEN
447*
448* Update Q := Q*Z**H
449*
450 CALL cunmr2( 'Right', 'Conjugate transpose', n, n, l, b,
451 $ ldb, tau, q, ldq, work, info )
452 END IF
453*
454* Clean up B
455*
456 CALL claset( 'Full', l, n-l, czero, czero, b, ldb )
457 DO 60 j = n - l + 1, n
458 DO 50 i = j - n + l + 1, l
459 b( i, j ) = czero
460 50 CONTINUE
461 60 CONTINUE
462*
463 END IF
464*
465* Let N-L L
466* A = ( A11 A12 ) M,
467*
468* then the following does the complete QR decomposition of A11:
469*
470* A11 = U*( 0 T12 )*P1**H
471* ( 0 0 )
472*
473 DO 70 i = 1, n - l
474 iwork( i ) = 0
475 70 CONTINUE
476 CALL cgeqp3( m, n-l, a, lda, iwork, tau, work, lwork, rwork,
477 $ info )
478*
479* Determine the effective rank of A11
480*
481 k = 0
482 DO 80 i = 1, min( m, n-l )
483 IF( abs( a( i, i ) ).GT.tola )
484 $ k = k + 1
485 80 CONTINUE
486*
487* Update A12 := U**H*A12, where A12 = A( 1:M, N-L+1:N )
488*
489 CALL cunm2r( 'Left', 'Conjugate transpose', m, l, min( m, n-l ),
490 $ a, lda, tau, a( 1, n-l+1 ), lda, work, info )
491*
492 IF( wantu ) THEN
493*
494* Copy the details of U, and form U
495*
496 CALL claset( 'Full', m, m, czero, czero, u, ldu )
497 IF( m.GT.1 )
498 $ CALL clacpy( 'Lower', m-1, n-l, a( 2, 1 ), lda, u( 2, 1 ),
499 $ ldu )
500 CALL cung2r( m, m, min( m, n-l ), u, ldu, tau, work, info )
501 END IF
502*
503 IF( wantq ) THEN
504*
505* Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1
506*
507 CALL clapmt( forwrd, n, n-l, q, ldq, iwork )
508 END IF
509*
510* Clean up A: set the strictly lower triangular part of
511* A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0.
512*
513 DO 100 j = 1, k - 1
514 DO 90 i = j + 1, k
515 a( i, j ) = czero
516 90 CONTINUE
517 100 CONTINUE
518 IF( m.GT.k )
519 $ CALL claset( 'Full', m-k, n-l, czero, czero, a( k+1, 1 ), lda )
520*
521 IF( n-l.GT.k ) THEN
522*
523* RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1
524*
525 CALL cgerq2( k, n-l, a, lda, tau, work, info )
526*
527 IF( wantq ) THEN
528*
529* Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**H
530*
531 CALL cunmr2( 'Right', 'Conjugate transpose', n, n-l, k, a,
532 $ lda, tau, q, ldq, work, info )
533 END IF
534*
535* Clean up A
536*
537 CALL claset( 'Full', k, n-l-k, czero, czero, a, lda )
538 DO 120 j = n - l - k + 1, n - l
539 DO 110 i = j - n + l + k + 1, k
540 a( i, j ) = czero
541 110 CONTINUE
542 120 CONTINUE
543*
544 END IF
545*
546 IF( m.GT.k ) THEN
547*
548* QR factorization of A( K+1:M,N-L+1:N )
549*
550 CALL cgeqr2( m-k, l, a( k+1, n-l+1 ), lda, tau, work, info )
551*
552 IF( wantu ) THEN
553*
554* Update U(:,K+1:M) := U(:,K+1:M)*U1
555*
556 CALL cunm2r( 'Right', 'No transpose', m, m-k, min( m-k, l ),
557 $ a( k+1, n-l+1 ), lda, tau, u( 1, k+1 ), ldu,
558 $ work, info )
559 END IF
560*
561* Clean up
562*
563 DO 140 j = n - l + 1, n
564 DO 130 i = j - n + k + l + 1, m
565 a( i, j ) = czero
566 130 CONTINUE
567 140 CONTINUE
568*
569 END IF
570*
571 work( 1 ) = cmplx( lwkopt )
572 RETURN
573*
574* End of CGGSVP3
575*
576 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cgeqp3(m, n, a, lda, jpvt, tau, work, lwork, rwork, info)
CGEQP3
Definition cgeqp3.f:159
subroutine cgeqr2(m, n, a, lda, tau, work, info)
CGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
Definition cgeqr2.f:130
subroutine cgerq2(m, n, a, lda, tau, work, info)
CGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm.
Definition cgerq2.f:123
subroutine cggsvp3(jobu, jobv, jobq, m, p, n, a, lda, b, ldb, tola, tolb, k, l, u, ldu, v, ldv, q, ldq, iwork, rwork, tau, work, lwork, info)
CGGSVP3
Definition cggsvp3.f:278
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
subroutine clapmt(forwrd, m, n, x, ldx, k)
CLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition clapmt.f:104
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
subroutine cung2r(m, n, k, a, lda, tau, work, info)
CUNG2R
Definition cung2r.f:114
subroutine cunm2r(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
CUNM2R multiplies a general matrix by the unitary matrix from a QR factorization determined by cgeqrf...
Definition cunm2r.f:159
subroutine cunmr2(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
CUNMR2 multiplies a general matrix by the unitary matrix from a RQ factorization determined by cgerqf...
Definition cunmr2.f:159