LAPACK  3.10.0
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 *> \ingroup complexOTHERsing
327 *
328 *> \par Contributors:
329 * ==================
330 *>
331 *> Ming Gu and Huan Ren, Computer Science Division, University of
332 *> California at Berkeley, USA
333 *>
334 * =====================================================================
335  SUBROUTINE cggsvd( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
336  $ LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK,
337  $ RWORK, IWORK, INFO )
338 *
339 * -- LAPACK driver routine --
340 * -- LAPACK is a software package provided by Univ. of Tennessee, --
341 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
342 *
343 * .. Scalar Arguments ..
344  CHARACTER JOBQ, JOBU, JOBV
345  INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
346 * ..
347 * .. Array Arguments ..
348  INTEGER IWORK( * )
349  REAL ALPHA( * ), BETA( * ), RWORK( * )
350  COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
351  $ u( ldu, * ), v( ldv, * ), work( * )
352 * ..
353 *
354 * =====================================================================
355 *
356 * .. Local Scalars ..
357  LOGICAL WANTQ, WANTU, WANTV
358  INTEGER I, IBND, ISUB, J, NCYCLE
359  REAL ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
360 * ..
361 * .. External Functions ..
362  LOGICAL LSAME
363  REAL CLANGE, SLAMCH
364  EXTERNAL lsame, clange, slamch
365 * ..
366 * .. External Subroutines ..
367  EXTERNAL cggsvp, ctgsja, scopy, xerbla
368 * ..
369 * .. Intrinsic Functions ..
370  INTRINSIC max, min
371 * ..
372 * .. Executable Statements ..
373 *
374 * Decode and test the input parameters
375 *
376  wantu = lsame( jobu, 'U' )
377  wantv = lsame( jobv, 'V' )
378  wantq = lsame( jobq, 'Q' )
379 *
380  info = 0
381  IF( .NOT.( wantu .OR. lsame( jobu, 'N' ) ) ) THEN
382  info = -1
383  ELSE IF( .NOT.( wantv .OR. lsame( jobv, 'N' ) ) ) THEN
384  info = -2
385  ELSE IF( .NOT.( wantq .OR. lsame( jobq, 'N' ) ) ) THEN
386  info = -3
387  ELSE IF( m.LT.0 ) THEN
388  info = -4
389  ELSE IF( n.LT.0 ) THEN
390  info = -5
391  ELSE IF( p.LT.0 ) THEN
392  info = -6
393  ELSE IF( lda.LT.max( 1, m ) ) THEN
394  info = -10
395  ELSE IF( ldb.LT.max( 1, p ) ) THEN
396  info = -12
397  ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) ) THEN
398  info = -16
399  ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) ) THEN
400  info = -18
401  ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
402  info = -20
403  END IF
404  IF( info.NE.0 ) THEN
405  CALL xerbla( 'CGGSVD', -info )
406  RETURN
407  END IF
408 *
409 * Compute the Frobenius norm of matrices A and B
410 *
411  anorm = clange( '1', m, n, a, lda, rwork )
412  bnorm = clange( '1', p, n, b, ldb, rwork )
413 *
414 * Get machine precision and set up threshold for determining
415 * the effective numerical rank of the matrices A and B.
416 *
417  ulp = slamch( 'Precision' )
418  unfl = slamch( 'Safe Minimum' )
419  tola = max( m, n )*max( anorm, unfl )*ulp
420  tolb = max( p, n )*max( bnorm, unfl )*ulp
421 *
422  CALL cggsvp( jobu, jobv, jobq, m, p, n, a, lda, b, ldb, tola,
423  $ tolb, k, l, u, ldu, v, ldv, q, ldq, iwork, rwork,
424  $ work, work( n+1 ), info )
425 *
426 * Compute the GSVD of two upper "triangular" matrices
427 *
428  CALL ctgsja( jobu, jobv, jobq, m, p, n, k, l, a, lda, b, ldb,
429  $ tola, tolb, alpha, beta, u, ldu, v, ldv, q, ldq,
430  $ work, ncycle, info )
431 *
432 * Sort the singular values and store the pivot indices in IWORK
433 * Copy ALPHA to RWORK, then sort ALPHA in RWORK
434 *
435  CALL scopy( n, alpha, 1, rwork, 1 )
436  ibnd = min( l, m-k )
437  DO 20 i = 1, ibnd
438 *
439 * Scan for largest ALPHA(K+I)
440 *
441  isub = i
442  smax = rwork( k+i )
443  DO 10 j = i + 1, ibnd
444  temp = rwork( k+j )
445  IF( temp.GT.smax ) THEN
446  isub = j
447  smax = temp
448  END IF
449  10 CONTINUE
450  IF( isub.NE.i ) THEN
451  rwork( k+isub ) = rwork( k+i )
452  rwork( k+i ) = smax
453  iwork( k+i ) = k + isub
454  ELSE
455  iwork( k+i ) = k + i
456  END IF
457  20 CONTINUE
458 *
459  RETURN
460 *
461 * End of CGGSVD
462 *
463  END
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:338
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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:262
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:379
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82