LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
cggsvd.f
Go to the documentation of this file.
1 *> \brief <b> CGGSVD computes the singular value decomposition (SVD) for OTHER matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CGGSVD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggsvd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggsvd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggsvd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
22 * LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK,
23 * RWORK, IWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBQ, JOBU, JOBV
27 * INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
28 * ..
29 * .. Array Arguments ..
30 * INTEGER IWORK( * )
31 * REAL ALPHA( * ), BETA( * ), RWORK( * )
32 * COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
33 * $ U( LDU, * ), V( LDV, * ), WORK( * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> This routine is deprecated and has been replaced by routine CGGSVD3.
43 *>
44 *> CGGSVD computes the generalized singular value decomposition (GSVD)
45 *> of an M-by-N complex matrix A and P-by-N complex matrix B:
46 *>
47 *> U**H*A*Q = D1*( 0 R ), V**H*B*Q = D2*( 0 R )
48 *>
49 *> where U, V and Q are unitary matrices.
50 *> Let K+L = the effective numerical rank of the
51 *> matrix (A**H,B**H)**H, then R is a (K+L)-by-(K+L) nonsingular upper
52 *> triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal"
53 *> matrices and of the following structures, respectively:
54 *>
55 *> If M-K-L >= 0,
56 *>
57 *> K L
58 *> D1 = K ( I 0 )
59 *> L ( 0 C )
60 *> M-K-L ( 0 0 )
61 *>
62 *> K L
63 *> D2 = L ( 0 S )
64 *> P-L ( 0 0 )
65 *>
66 *> N-K-L K L
67 *> ( 0 R ) = K ( 0 R11 R12 )
68 *> L ( 0 0 R22 )
69 *>
70 *> where
71 *>
72 *> C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
73 *> S = diag( BETA(K+1), ... , BETA(K+L) ),
74 *> C**2 + S**2 = I.
75 *>
76 *> R is stored in A(1:K+L,N-K-L+1:N) on exit.
77 *>
78 *> If M-K-L < 0,
79 *>
80 *> K M-K K+L-M
81 *> D1 = K ( I 0 0 )
82 *> M-K ( 0 C 0 )
83 *>
84 *> K M-K K+L-M
85 *> D2 = M-K ( 0 S 0 )
86 *> K+L-M ( 0 0 I )
87 *> P-L ( 0 0 0 )
88 *>
89 *> N-K-L K M-K K+L-M
90 *> ( 0 R ) = K ( 0 R11 R12 R13 )
91 *> M-K ( 0 0 R22 R23 )
92 *> K+L-M ( 0 0 0 R33 )
93 *>
94 *> where
95 *>
96 *> C = diag( ALPHA(K+1), ... , ALPHA(M) ),
97 *> S = diag( BETA(K+1), ... , BETA(M) ),
98 *> C**2 + S**2 = I.
99 *>
100 *> (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
101 *> ( 0 R22 R23 )
102 *> in B(M-K+1:L,N+M-K-L+1:N) on exit.
103 *>
104 *> The routine computes C, S, R, and optionally the unitary
105 *> transformation matrices U, V and Q.
106 *>
107 *> In particular, if B is an N-by-N nonsingular matrix, then the GSVD of
108 *> A and B implicitly gives the SVD of A*inv(B):
109 *> A*inv(B) = U*(D1*inv(D2))*V**H.
110 *> If ( A**H,B**H)**H has orthnormal columns, then the GSVD of A and B is also
111 *> equal to the CS decomposition of A and B. Furthermore, the GSVD can
112 *> be used to derive the solution of the eigenvalue problem:
113 *> A**H*A x = lambda* B**H*B x.
114 *> In some literature, the GSVD of A and B is presented in the form
115 *> U**H*A*X = ( 0 D1 ), V**H*B*X = ( 0 D2 )
116 *> where U and V are orthogonal and X is nonsingular, and D1 and D2 are
117 *> ``diagonal''. The former GSVD form can be converted to the latter
118 *> form by taking the nonsingular matrix X as
119 *>
120 *> X = Q*( I 0 )
121 *> ( 0 inv(R) )
122 *> \endverbatim
123 *
124 * Arguments:
125 * ==========
126 *
127 *> \param[in] JOBU
128 *> \verbatim
129 *> JOBU is CHARACTER*1
130 *> = 'U': Unitary matrix U is computed;
131 *> = 'N': U is not computed.
132 *> \endverbatim
133 *>
134 *> \param[in] JOBV
135 *> \verbatim
136 *> JOBV is CHARACTER*1
137 *> = 'V': Unitary matrix V is computed;
138 *> = 'N': V is not computed.
139 *> \endverbatim
140 *>
141 *> \param[in] JOBQ
142 *> \verbatim
143 *> JOBQ is CHARACTER*1
144 *> = 'Q': Unitary matrix Q is computed;
145 *> = 'N': Q is not computed.
146 *> \endverbatim
147 *>
148 *> \param[in] M
149 *> \verbatim
150 *> M is INTEGER
151 *> The number of rows of the matrix A. M >= 0.
152 *> \endverbatim
153 *>
154 *> \param[in] N
155 *> \verbatim
156 *> N is INTEGER
157 *> The number of columns of the matrices A and B. N >= 0.
158 *> \endverbatim
159 *>
160 *> \param[in] P
161 *> \verbatim
162 *> P is INTEGER
163 *> The number of rows of the matrix B. P >= 0.
164 *> \endverbatim
165 *>
166 *> \param[out] K
167 *> \verbatim
168 *> K is INTEGER
169 *> \endverbatim
170 *>
171 *> \param[out] L
172 *> \verbatim
173 *> L is INTEGER
174 *>
175 *> On exit, K and L specify the dimension of the subblocks
176 *> described in Purpose.
177 *> K + L = effective numerical rank of (A**H,B**H)**H.
178 *> \endverbatim
179 *>
180 *> \param[in,out] A
181 *> \verbatim
182 *> A is COMPLEX array, dimension (LDA,N)
183 *> On entry, the M-by-N matrix A.
184 *> On exit, A contains the triangular matrix R, or part of R.
185 *> See Purpose for details.
186 *> \endverbatim
187 *>
188 *> \param[in] LDA
189 *> \verbatim
190 *> LDA is INTEGER
191 *> The leading dimension of the array A. LDA >= max(1,M).
192 *> \endverbatim
193 *>
194 *> \param[in,out] B
195 *> \verbatim
196 *> B is COMPLEX array, dimension (LDB,N)
197 *> On entry, the P-by-N matrix B.
198 *> On exit, B contains part of the triangular matrix R if
199 *> M-K-L < 0. See Purpose for details.
200 *> \endverbatim
201 *>
202 *> \param[in] LDB
203 *> \verbatim
204 *> LDB is INTEGER
205 *> The leading dimension of the array B. LDB >= max(1,P).
206 *> \endverbatim
207 *>
208 *> \param[out] ALPHA
209 *> \verbatim
210 *> ALPHA is REAL array, dimension (N)
211 *> \endverbatim
212 *>
213 *> \param[out] BETA
214 *> \verbatim
215 *> BETA is REAL array, dimension (N)
216 *>
217 *> On exit, ALPHA and BETA contain the generalized singular
218 *> value pairs of A and B;
219 *> ALPHA(1:K) = 1,
220 *> BETA(1:K) = 0,
221 *> and if M-K-L >= 0,
222 *> ALPHA(K+1:K+L) = C,
223 *> BETA(K+1:K+L) = S,
224 *> or if M-K-L < 0,
225 *> ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0
226 *> BETA(K+1:M) =S, BETA(M+1:K+L) =1
227 *> and
228 *> ALPHA(K+L+1:N) = 0
229 *> BETA(K+L+1:N) = 0
230 *> \endverbatim
231 *>
232 *> \param[out] U
233 *> \verbatim
234 *> U is COMPLEX array, dimension (LDU,M)
235 *> If JOBU = 'U', U contains the M-by-M unitary matrix U.
236 *> If JOBU = 'N', U is not referenced.
237 *> \endverbatim
238 *>
239 *> \param[in] LDU
240 *> \verbatim
241 *> LDU is INTEGER
242 *> The leading dimension of the array U. LDU >= max(1,M) if
243 *> JOBU = 'U'; LDU >= 1 otherwise.
244 *> \endverbatim
245 *>
246 *> \param[out] V
247 *> \verbatim
248 *> V is COMPLEX array, dimension (LDV,P)
249 *> If JOBV = 'V', V contains the P-by-P unitary matrix V.
250 *> If JOBV = 'N', V is not referenced.
251 *> \endverbatim
252 *>
253 *> \param[in] LDV
254 *> \verbatim
255 *> LDV is INTEGER
256 *> The leading dimension of the array V. LDV >= max(1,P) if
257 *> JOBV = 'V'; LDV >= 1 otherwise.
258 *> \endverbatim
259 *>
260 *> \param[out] Q
261 *> \verbatim
262 *> Q is COMPLEX array, dimension (LDQ,N)
263 *> If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q.
264 *> If JOBQ = 'N', Q is not referenced.
265 *> \endverbatim
266 *>
267 *> \param[in] LDQ
268 *> \verbatim
269 *> LDQ is INTEGER
270 *> The leading dimension of the array Q. LDQ >= max(1,N) if
271 *> JOBQ = 'Q'; LDQ >= 1 otherwise.
272 *> \endverbatim
273 *>
274 *> \param[out] WORK
275 *> \verbatim
276 *> WORK is COMPLEX array, dimension (max(3*N,M,P)+N)
277 *> \endverbatim
278 *>
279 *> \param[out] RWORK
280 *> \verbatim
281 *> RWORK is REAL array, dimension (2*N)
282 *> \endverbatim
283 *>
284 *> \param[out] IWORK
285 *> \verbatim
286 *> IWORK is INTEGER array, dimension (N)
287 *> On exit, IWORK stores the sorting information. More
288 *> precisely, the following loop will sort ALPHA
289 *> for I = K+1, min(M,K+L)
290 *> swap ALPHA(I) and ALPHA(IWORK(I))
291 *> endfor
292 *> such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N).
293 *> \endverbatim
294 *>
295 *> \param[out] INFO
296 *> \verbatim
297 *> INFO is INTEGER
298 *> = 0: successful exit.
299 *> < 0: if INFO = -i, the i-th argument had an illegal value.
300 *> > 0: if INFO = 1, the Jacobi-type procedure failed to
301 *> converge. For further details, see subroutine CTGSJA.
302 *> \endverbatim
303 *
304 *> \par Internal Parameters:
305 * =========================
306 *>
307 *> \verbatim
308 *> TOLA REAL
309 *> TOLB REAL
310 *> TOLA and TOLB are the thresholds to determine the effective
311 *> rank of (A**H,B**H)**H. Generally, they are set to
312 *> TOLA = MAX(M,N)*norm(A)*MACHEPS,
313 *> TOLB = MAX(P,N)*norm(B)*MACHEPS.
314 *> The size of TOLA and TOLB may affect the size of backward
315 *> errors of the decomposition.
316 *> \endverbatim
317 *
318 * Authors:
319 * ========
320 *
321 *> \author Univ. of Tennessee
322 *> \author Univ. of California Berkeley
323 *> \author Univ. of Colorado Denver
324 *> \author NAG Ltd.
325 *
326 *> \date November 2011
327 *
328 *> \ingroup complexOTHERsing
329 *
330 *> \par Contributors:
331 * ==================
332 *>
333 *> Ming Gu and Huan Ren, Computer Science Division, University of
334 *> California at Berkeley, USA
335 *>
336 * =====================================================================
337  SUBROUTINE cggsvd( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
338  $ ldb, alpha, beta, u, ldu, v, ldv, q, ldq, work,
339  $ rwork, iwork, info )
340 *
341 * -- LAPACK driver routine (version 3.4.0) --
342 * -- LAPACK is a software package provided by Univ. of Tennessee, --
343 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
344 * November 2011
345 *
346 * .. Scalar Arguments ..
347  CHARACTER JOBQ, JOBU, JOBV
348  INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
349 * ..
350 * .. Array Arguments ..
351  INTEGER IWORK( * )
352  REAL ALPHA( * ), BETA( * ), RWORK( * )
353  COMPLEX A( lda, * ), B( ldb, * ), Q( ldq, * ),
354  $ u( ldu, * ), v( ldv, * ), work( * )
355 * ..
356 *
357 * =====================================================================
358 *
359 * .. Local Scalars ..
360  LOGICAL WANTQ, WANTU, WANTV
361  INTEGER I, IBND, ISUB, J, NCYCLE
362  REAL ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
363 * ..
364 * .. External Functions ..
365  LOGICAL LSAME
366  REAL CLANGE, SLAMCH
367  EXTERNAL lsame, clange, slamch
368 * ..
369 * .. External Subroutines ..
370  EXTERNAL cggsvp, ctgsja, scopy, xerbla
371 * ..
372 * .. Intrinsic Functions ..
373  INTRINSIC max, min
374 * ..
375 * .. Executable Statements ..
376 *
377 * Decode and test the input parameters
378 *
379  wantu = lsame( jobu, 'U' )
380  wantv = lsame( jobv, 'V' )
381  wantq = lsame( jobq, 'Q' )
382 *
383  info = 0
384  IF( .NOT.( wantu .OR. lsame( jobu, 'N' ) ) ) THEN
385  info = -1
386  ELSE IF( .NOT.( wantv .OR. lsame( jobv, 'N' ) ) ) THEN
387  info = -2
388  ELSE IF( .NOT.( wantq .OR. lsame( jobq, 'N' ) ) ) THEN
389  info = -3
390  ELSE IF( m.LT.0 ) THEN
391  info = -4
392  ELSE IF( n.LT.0 ) THEN
393  info = -5
394  ELSE IF( p.LT.0 ) THEN
395  info = -6
396  ELSE IF( lda.LT.max( 1, m ) ) THEN
397  info = -10
398  ELSE IF( ldb.LT.max( 1, p ) ) THEN
399  info = -12
400  ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) ) THEN
401  info = -16
402  ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) ) THEN
403  info = -18
404  ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
405  info = -20
406  END IF
407  IF( info.NE.0 ) THEN
408  CALL xerbla( 'CGGSVD', -info )
409  RETURN
410  END IF
411 *
412 * Compute the Frobenius norm of matrices A and B
413 *
414  anorm = clange( '1', m, n, a, lda, rwork )
415  bnorm = clange( '1', p, n, b, ldb, rwork )
416 *
417 * Get machine precision and set up threshold for determining
418 * the effective numerical rank of the matrices A and B.
419 *
420  ulp = slamch( 'Precision' )
421  unfl = slamch( 'Safe Minimum' )
422  tola = max( m, n )*max( anorm, unfl )*ulp
423  tolb = max( p, n )*max( bnorm, unfl )*ulp
424 *
425  CALL cggsvp( jobu, jobv, jobq, m, p, n, a, lda, b, ldb, tola,
426  $ tolb, k, l, u, ldu, v, ldv, q, ldq, iwork, rwork,
427  $ work, work( n+1 ), info )
428 *
429 * Compute the GSVD of two upper "triangular" matrices
430 *
431  CALL ctgsja( jobu, jobv, jobq, m, p, n, k, l, a, lda, b, ldb,
432  $ tola, tolb, alpha, beta, u, ldu, v, ldv, q, ldq,
433  $ work, ncycle, info )
434 *
435 * Sort the singular values and store the pivot indices in IWORK
436 * Copy ALPHA to RWORK, then sort ALPHA in RWORK
437 *
438  CALL scopy( n, alpha, 1, rwork, 1 )
439  ibnd = min( l, m-k )
440  DO 20 i = 1, ibnd
441 *
442 * Scan for largest ALPHA(K+I)
443 *
444  isub = i
445  smax = rwork( k+i )
446  DO 10 j = i + 1, ibnd
447  temp = rwork( k+j )
448  IF( temp.GT.smax ) THEN
449  isub = j
450  smax = temp
451  END IF
452  10 CONTINUE
453  IF( isub.NE.i ) THEN
454  rwork( k+isub ) = rwork( k+i )
455  rwork( k+i ) = smax
456  iwork( k+i ) = k + isub
457  ELSE
458  iwork( k+i ) = k + i
459  END IF
460  20 CONTINUE
461 *
462  RETURN
463 *
464 * End of CGGSVD
465 *
466  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cggsvd(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, RWORK, IWORK, INFO)
CGGSVD computes the singular value decomposition (SVD) for OTHER matrices
Definition: cggsvd.f:340
subroutine cggsvp(JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK, TAU, WORK, INFO)
CGGSVP
Definition: cggsvp.f:264
subroutine ctgsja(JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, NCYCLE, INFO)
CTGSJA
Definition: ctgsja.f:381
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53