LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
zggsvd.f
Go to the documentation of this file.
1 *> \brief <b> ZGGSVD 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 ZGGSVD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggsvd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggsvd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggsvd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZGGSVD( 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 * DOUBLE PRECISION ALPHA( * ), BETA( * ), RWORK( * )
32 * COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
33 * $ U( LDU, * ), V( LDV, * ), WORK( * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> ZGGSVD computes the generalized singular value decomposition (GSVD)
43 *> of an M-by-N complex matrix A and P-by-N complex matrix B:
44 *>
45 *> U**H*A*Q = D1*( 0 R ), V**H*B*Q = D2*( 0 R )
46 *>
47 *> where U, V and Q are unitary matrices.
48 *> Let K+L = the effective numerical rank of the
49 *> matrix (A**H,B**H)**H, then R is a (K+L)-by-(K+L) nonsingular upper
50 *> triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal"
51 *> matrices and of the following structures, respectively:
52 *>
53 *> If M-K-L >= 0,
54 *>
55 *> K L
56 *> D1 = K ( I 0 )
57 *> L ( 0 C )
58 *> M-K-L ( 0 0 )
59 *>
60 *> K L
61 *> D2 = L ( 0 S )
62 *> P-L ( 0 0 )
63 *>
64 *> N-K-L K L
65 *> ( 0 R ) = K ( 0 R11 R12 )
66 *> L ( 0 0 R22 )
67 *> where
68 *>
69 *> C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
70 *> S = diag( BETA(K+1), ... , BETA(K+L) ),
71 *> C**2 + S**2 = I.
72 *>
73 *> R is stored in A(1:K+L,N-K-L+1:N) on exit.
74 *>
75 *> If M-K-L < 0,
76 *>
77 *> K M-K K+L-M
78 *> D1 = K ( I 0 0 )
79 *> M-K ( 0 C 0 )
80 *>
81 *> K M-K K+L-M
82 *> D2 = M-K ( 0 S 0 )
83 *> K+L-M ( 0 0 I )
84 *> P-L ( 0 0 0 )
85 *>
86 *> N-K-L K M-K K+L-M
87 *> ( 0 R ) = K ( 0 R11 R12 R13 )
88 *> M-K ( 0 0 R22 R23 )
89 *> K+L-M ( 0 0 0 R33 )
90 *>
91 *> where
92 *>
93 *> C = diag( ALPHA(K+1), ... , ALPHA(M) ),
94 *> S = diag( BETA(K+1), ... , BETA(M) ),
95 *> C**2 + S**2 = I.
96 *>
97 *> (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
98 *> ( 0 R22 R23 )
99 *> in B(M-K+1:L,N+M-K-L+1:N) on exit.
100 *>
101 *> The routine computes C, S, R, and optionally the unitary
102 *> transformation matrices U, V and Q.
103 *>
104 *> In particular, if B is an N-by-N nonsingular matrix, then the GSVD of
105 *> A and B implicitly gives the SVD of A*inv(B):
106 *> A*inv(B) = U*(D1*inv(D2))*V**H.
107 *> If ( A**H,B**H)**H has orthnormal columns, then the GSVD of A and B is also
108 *> equal to the CS decomposition of A and B. Furthermore, the GSVD can
109 *> be used to derive the solution of the eigenvalue problem:
110 *> A**H*A x = lambda* B**H*B x.
111 *> In some literature, the GSVD of A and B is presented in the form
112 *> U**H*A*X = ( 0 D1 ), V**H*B*X = ( 0 D2 )
113 *> where U and V are orthogonal and X is nonsingular, and D1 and D2 are
114 *> ``diagonal''. The former GSVD form can be converted to the latter
115 *> form by taking the nonsingular matrix X as
116 *>
117 *> X = Q*( I 0 )
118 *> ( 0 inv(R) )
119 *> \endverbatim
120 *
121 * Arguments:
122 * ==========
123 *
124 *> \param[in] JOBU
125 *> \verbatim
126 *> JOBU is CHARACTER*1
127 *> = 'U': Unitary matrix U is computed;
128 *> = 'N': U is not computed.
129 *> \endverbatim
130 *>
131 *> \param[in] JOBV
132 *> \verbatim
133 *> JOBV is CHARACTER*1
134 *> = 'V': Unitary matrix V is computed;
135 *> = 'N': V is not computed.
136 *> \endverbatim
137 *>
138 *> \param[in] JOBQ
139 *> \verbatim
140 *> JOBQ is CHARACTER*1
141 *> = 'Q': Unitary matrix Q is computed;
142 *> = 'N': Q is not computed.
143 *> \endverbatim
144 *>
145 *> \param[in] M
146 *> \verbatim
147 *> M is INTEGER
148 *> The number of rows of the matrix A. M >= 0.
149 *> \endverbatim
150 *>
151 *> \param[in] N
152 *> \verbatim
153 *> N is INTEGER
154 *> The number of columns of the matrices A and B. N >= 0.
155 *> \endverbatim
156 *>
157 *> \param[in] P
158 *> \verbatim
159 *> P is INTEGER
160 *> The number of rows of the matrix B. P >= 0.
161 *> \endverbatim
162 *>
163 *> \param[out] K
164 *> \verbatim
165 *> K is INTEGER
166 *> \endverbatim
167 *>
168 *> \param[out] L
169 *> \verbatim
170 *> L is INTEGER
171 *>
172 *> On exit, K and L specify the dimension of the subblocks
173 *> described in Purpose.
174 *> K + L = effective numerical rank of (A**H,B**H)**H.
175 *> \endverbatim
176 *>
177 *> \param[in,out] A
178 *> \verbatim
179 *> A is COMPLEX*16 array, dimension (LDA,N)
180 *> On entry, the M-by-N matrix A.
181 *> On exit, A contains the triangular matrix R, or part of R.
182 *> See Purpose for details.
183 *> \endverbatim
184 *>
185 *> \param[in] LDA
186 *> \verbatim
187 *> LDA is INTEGER
188 *> The leading dimension of the array A. LDA >= max(1,M).
189 *> \endverbatim
190 *>
191 *> \param[in,out] B
192 *> \verbatim
193 *> B is COMPLEX*16 array, dimension (LDB,N)
194 *> On entry, the P-by-N matrix B.
195 *> On exit, B contains part of the triangular matrix R if
196 *> M-K-L < 0. See Purpose for details.
197 *> \endverbatim
198 *>
199 *> \param[in] LDB
200 *> \verbatim
201 *> LDB is INTEGER
202 *> The leading dimension of the array B. LDB >= max(1,P).
203 *> \endverbatim
204 *>
205 *> \param[out] ALPHA
206 *> \verbatim
207 *> ALPHA is DOUBLE PRECISION array, dimension (N)
208 *> \endverbatim
209 *>
210 *> \param[out] BETA
211 *> \verbatim
212 *> BETA is DOUBLE PRECISION array, dimension (N)
213 *>
214 *> On exit, ALPHA and BETA contain the generalized singular
215 *> value pairs of A and B;
216 *> ALPHA(1:K) = 1,
217 *> BETA(1:K) = 0,
218 *> and if M-K-L >= 0,
219 *> ALPHA(K+1:K+L) = C,
220 *> BETA(K+1:K+L) = S,
221 *> or if M-K-L < 0,
222 *> ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0
223 *> BETA(K+1:M) =S, BETA(M+1:K+L) =1
224 *> and
225 *> ALPHA(K+L+1:N) = 0
226 *> BETA(K+L+1:N) = 0
227 *> \endverbatim
228 *>
229 *> \param[out] U
230 *> \verbatim
231 *> U is COMPLEX*16 array, dimension (LDU,M)
232 *> If JOBU = 'U', U contains the M-by-M unitary matrix U.
233 *> If JOBU = 'N', U is not referenced.
234 *> \endverbatim
235 *>
236 *> \param[in] LDU
237 *> \verbatim
238 *> LDU is INTEGER
239 *> The leading dimension of the array U. LDU >= max(1,M) if
240 *> JOBU = 'U'; LDU >= 1 otherwise.
241 *> \endverbatim
242 *>
243 *> \param[out] V
244 *> \verbatim
245 *> V is COMPLEX*16 array, dimension (LDV,P)
246 *> If JOBV = 'V', V contains the P-by-P unitary matrix V.
247 *> If JOBV = 'N', V is not referenced.
248 *> \endverbatim
249 *>
250 *> \param[in] LDV
251 *> \verbatim
252 *> LDV is INTEGER
253 *> The leading dimension of the array V. LDV >= max(1,P) if
254 *> JOBV = 'V'; LDV >= 1 otherwise.
255 *> \endverbatim
256 *>
257 *> \param[out] Q
258 *> \verbatim
259 *> Q is COMPLEX*16 array, dimension (LDQ,N)
260 *> If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q.
261 *> If JOBQ = 'N', Q is not referenced.
262 *> \endverbatim
263 *>
264 *> \param[in] LDQ
265 *> \verbatim
266 *> LDQ is INTEGER
267 *> The leading dimension of the array Q. LDQ >= max(1,N) if
268 *> JOBQ = 'Q'; LDQ >= 1 otherwise.
269 *> \endverbatim
270 *>
271 *> \param[out] WORK
272 *> \verbatim
273 *> WORK is COMPLEX*16 array, dimension (max(3*N,M,P)+N)
274 *> \endverbatim
275 *>
276 *> \param[out] RWORK
277 *> \verbatim
278 *> RWORK is DOUBLE PRECISION array, dimension (2*N)
279 *> \endverbatim
280 *>
281 *> \param[out] IWORK
282 *> \verbatim
283 *> IWORK is INTEGER array, dimension (N)
284 *> On exit, IWORK stores the sorting information. More
285 *> precisely, the following loop will sort ALPHA
286 *> for I = K+1, min(M,K+L)
287 *> swap ALPHA(I) and ALPHA(IWORK(I))
288 *> endfor
289 *> such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N).
290 *> \endverbatim
291 *>
292 *> \param[out] INFO
293 *> \verbatim
294 *> INFO is INTEGER
295 *> = 0: successful exit.
296 *> < 0: if INFO = -i, the i-th argument had an illegal value.
297 *> > 0: if INFO = 1, the Jacobi-type procedure failed to
298 *> converge. For further details, see subroutine ZTGSJA.
299 *> \endverbatim
300 *
301 *> \par Internal Parameters:
302 * =========================
303 *>
304 *> \verbatim
305 *> TOLA DOUBLE PRECISION
306 *> TOLB DOUBLE PRECISION
307 *> TOLA and TOLB are the thresholds to determine the effective
308 *> rank of (A**H,B**H)**H. Generally, they are set to
309 *> TOLA = MAX(M,N)*norm(A)*MAZHEPS,
310 *> TOLB = MAX(P,N)*norm(B)*MAZHEPS.
311 *> The size of TOLA and TOLB may affect the size of backward
312 *> errors of the decomposition.
313 *> \endverbatim
314 *
315 * Authors:
316 * ========
317 *
318 *> \author Univ. of Tennessee
319 *> \author Univ. of California Berkeley
320 *> \author Univ. of Colorado Denver
321 *> \author NAG Ltd.
322 *
323 *> \date November 2011
324 *
325 *> \ingroup complex16OTHERsing
326 *
327 *> \par Contributors:
328 * ==================
329 *>
330 *> Ming Gu and Huan Ren, Computer Science Division, University of
331 *> California at Berkeley, USA
332 *>
333 * =====================================================================
334  SUBROUTINE zggsvd( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
335  $ ldb, alpha, beta, u, ldu, v, ldv, q, ldq, work,
336  $ rwork, iwork, info )
337 *
338 * -- LAPACK driver routine (version 3.4.0) --
339 * -- LAPACK is a software package provided by Univ. of Tennessee, --
340 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
341 * November 2011
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  DOUBLE PRECISION ALPHA( * ), BETA( * ), RWORK( * )
350  COMPLEX*16 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  DOUBLE PRECISION ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
360 * ..
361 * .. External Functions ..
362  LOGICAL LSAME
363  DOUBLE PRECISION DLAMCH, ZLANGE
364  EXTERNAL lsame, dlamch, zlange
365 * ..
366 * .. External Subroutines ..
367  EXTERNAL dcopy, xerbla, zggsvp, ztgsja
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( 'ZGGSVD', -info )
406  RETURN
407  END IF
408 *
409 * Compute the Frobenius norm of matrices A and B
410 *
411  anorm = zlange( '1', m, n, a, lda, rwork )
412  bnorm = zlange( '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 = dlamch( 'Precision' )
418  unfl = dlamch( 'Safe Minimum' )
419  tola = max( m, n )*max( anorm, unfl )*ulp
420  tolb = max( p, n )*max( bnorm, unfl )*ulp
421 *
422  CALL zggsvp( 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 ztgsja( 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 dcopy( 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 ZGGSVD
462 *
463  END
subroutine zggsvp(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)
ZGGSVP
Definition: zggsvp.f:265
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zggsvd(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, RWORK, IWORK, INFO)
ZGGSVD computes the singular value decomposition (SVD) for OTHER matrices
Definition: zggsvd.f:337
subroutine ztgsja(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)
ZTGSJA
Definition: ztgsja.f:381