LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zgsvts3.f
Go to the documentation of this file.
1*> \brief \b ZGSVTS3
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE ZGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
12* LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK,
13* LWORK, RWORK, RESULT )
14*
15* .. Scalar Arguments ..
16* INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
17* ..
18* .. Array Arguments ..
19* INTEGER IWORK( * )
20* DOUBLE PRECISION ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
21* COMPLEX*16 A( LDA, * ), AF( LDA, * ), B( LDB, * ),
22* $ BF( LDB, * ), Q( LDQ, * ), R( LDR, * ),
23* $ U( LDU, * ), V( LDV, * ), WORK( LWORK )
24* ..
25*
26*
27*> \par Purpose:
28* =============
29*>
30*> \verbatim
31*>
32*> ZGSVTS3 tests ZGGSVD3, which computes the GSVD of an M-by-N matrix A
33*> and a P-by-N matrix B:
34*> U'*A*Q = D1*R and V'*B*Q = D2*R.
35*> \endverbatim
36*
37* Arguments:
38* ==========
39*
40*> \param[in] M
41*> \verbatim
42*> M is INTEGER
43*> The number of rows of the matrix A. M >= 0.
44*> \endverbatim
45*>
46*> \param[in] P
47*> \verbatim
48*> P is INTEGER
49*> The number of rows of the matrix B. P >= 0.
50*> \endverbatim
51*>
52*> \param[in] N
53*> \verbatim
54*> N is INTEGER
55*> The number of columns of the matrices A and B. N >= 0.
56*> \endverbatim
57*>
58*> \param[in] A
59*> \verbatim
60*> A is COMPLEX*16 array, dimension (LDA,M)
61*> The M-by-N matrix A.
62*> \endverbatim
63*>
64*> \param[out] AF
65*> \verbatim
66*> AF is COMPLEX*16 array, dimension (LDA,N)
67*> Details of the GSVD of A and B, as returned by ZGGSVD3,
68*> see ZGGSVD3 for further details.
69*> \endverbatim
70*>
71*> \param[in] LDA
72*> \verbatim
73*> LDA is INTEGER
74*> The leading dimension of the arrays A and AF.
75*> LDA >= max( 1,M ).
76*> \endverbatim
77*>
78*> \param[in] B
79*> \verbatim
80*> B is COMPLEX*16 array, dimension (LDB,P)
81*> On entry, the P-by-N matrix B.
82*> \endverbatim
83*>
84*> \param[out] BF
85*> \verbatim
86*> BF is COMPLEX*16 array, dimension (LDB,N)
87*> Details of the GSVD of A and B, as returned by ZGGSVD3,
88*> see ZGGSVD3 for further details.
89*> \endverbatim
90*>
91*> \param[in] LDB
92*> \verbatim
93*> LDB is INTEGER
94*> The leading dimension of the arrays B and BF.
95*> LDB >= max(1,P).
96*> \endverbatim
97*>
98*> \param[out] U
99*> \verbatim
100*> U is COMPLEX*16 array, dimension(LDU,M)
101*> The M by M unitary matrix U.
102*> \endverbatim
103*>
104*> \param[in] LDU
105*> \verbatim
106*> LDU is INTEGER
107*> The leading dimension of the array U. LDU >= max(1,M).
108*> \endverbatim
109*>
110*> \param[out] V
111*> \verbatim
112*> V is COMPLEX*16 array, dimension(LDV,M)
113*> The P by P unitary matrix V.
114*> \endverbatim
115*>
116*> \param[in] LDV
117*> \verbatim
118*> LDV is INTEGER
119*> The leading dimension of the array V. LDV >= max(1,P).
120*> \endverbatim
121*>
122*> \param[out] Q
123*> \verbatim
124*> Q is COMPLEX*16 array, dimension(LDQ,N)
125*> The N by N unitary matrix Q.
126*> \endverbatim
127*>
128*> \param[in] LDQ
129*> \verbatim
130*> LDQ is INTEGER
131*> The leading dimension of the array Q. LDQ >= max(1,N).
132*> \endverbatim
133*>
134*> \param[out] ALPHA
135*> \verbatim
136*> ALPHA is DOUBLE PRECISION array, dimension (N)
137*> \endverbatim
138*>
139*> \param[out] BETA
140*> \verbatim
141*> BETA is DOUBLE PRECISION array, dimension (N)
142*>
143*> The generalized singular value pairs of A and B, the
144*> ``diagonal'' matrices D1 and D2 are constructed from
145*> ALPHA and BETA, see subroutine ZGGSVD3 for details.
146*> \endverbatim
147*>
148*> \param[out] R
149*> \verbatim
150*> R is COMPLEX*16 array, dimension(LDQ,N)
151*> The upper triangular matrix R.
152*> \endverbatim
153*>
154*> \param[in] LDR
155*> \verbatim
156*> LDR is INTEGER
157*> The leading dimension of the array R. LDR >= max(1,N).
158*> \endverbatim
159*>
160*> \param[out] IWORK
161*> \verbatim
162*> IWORK is INTEGER array, dimension (N)
163*> \endverbatim
164*>
165*> \param[out] WORK
166*> \verbatim
167*> WORK is COMPLEX*16 array, dimension (LWORK)
168*> \endverbatim
169*>
170*> \param[in] LWORK
171*> \verbatim
172*> LWORK is INTEGER
173*> The dimension of the array WORK,
174*> LWORK >= max(M,P,N)*max(M,P,N).
175*> \endverbatim
176*>
177*> \param[out] RWORK
178*> \verbatim
179*> RWORK is DOUBLE PRECISION array, dimension (max(M,P,N))
180*> \endverbatim
181*>
182*> \param[out] RESULT
183*> \verbatim
184*> RESULT is DOUBLE PRECISION array, dimension (6)
185*> The test ratios:
186*> RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP)
187*> RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP)
188*> RESULT(3) = norm( I - U'*U ) / ( M*ULP )
189*> RESULT(4) = norm( I - V'*V ) / ( P*ULP )
190*> RESULT(5) = norm( I - Q'*Q ) / ( N*ULP )
191*> RESULT(6) = 0 if ALPHA is in decreasing order;
192*> = ULPINV otherwise.
193*> \endverbatim
194*
195* Authors:
196* ========
197*
198*> \author Univ. of Tennessee
199*> \author Univ. of California Berkeley
200*> \author Univ. of Colorado Denver
201*> \author NAG Ltd.
202*
203*> \ingroup complex16_eig
204*
205* =====================================================================
206 SUBROUTINE zgsvts3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
207 $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK,
208 $ LWORK, RWORK, RESULT )
209*
210* -- LAPACK test routine --
211* -- LAPACK is a software package provided by Univ. of Tennessee, --
212* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213*
214* .. Scalar Arguments ..
215 INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
216* ..
217* .. Array Arguments ..
218 INTEGER IWORK( * )
219 DOUBLE PRECISION ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
220 COMPLEX*16 A( LDA, * ), AF( LDA, * ), B( LDB, * ),
221 $ bf( ldb, * ), q( ldq, * ), r( ldr, * ),
222 $ u( ldu, * ), v( ldv, * ), work( lwork )
223* ..
224*
225* =====================================================================
226*
227* .. Parameters ..
228 DOUBLE PRECISION ZERO, ONE
229 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
230 COMPLEX*16 CZERO, CONE
231 parameter( czero = ( 0.0d+0, 0.0d+0 ),
232 $ cone = ( 1.0d+0, 0.0d+0 ) )
233* ..
234* .. Local Scalars ..
235 INTEGER I, INFO, J, K, L
236 DOUBLE PRECISION ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
237* ..
238* .. External Functions ..
239 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE
240 EXTERNAL DLAMCH, ZLANGE, ZLANHE
241* ..
242* .. External Subroutines ..
243 EXTERNAL dcopy, zgemm, zggsvd3, zherk, zlacpy, zlaset
244* ..
245* .. Intrinsic Functions ..
246 INTRINSIC dble, max, min
247* ..
248* .. Executable Statements ..
249*
250 ulp = dlamch( 'Precision' )
251 ulpinv = one / ulp
252 unfl = dlamch( 'Safe minimum' )
253*
254* Copy the matrix A to the array AF.
255*
256 CALL zlacpy( 'Full', m, n, a, lda, af, lda )
257 CALL zlacpy( 'Full', p, n, b, ldb, bf, ldb )
258*
259 anorm = max( zlange( '1', m, n, a, lda, rwork ), unfl )
260 bnorm = max( zlange( '1', p, n, b, ldb, rwork ), unfl )
261*
262* Factorize the matrices A and B in the arrays AF and BF.
263*
264 CALL zggsvd3( 'U', 'V', 'Q', m, n, p, k, l, af, lda, bf, ldb,
265 $ alpha, beta, u, ldu, v, ldv, q, ldq, work, lwork,
266 $ rwork, iwork, info )
267*
268* Copy R
269*
270 DO 20 i = 1, min( k+l, m )
271 DO 10 j = i, k + l
272 r( i, j ) = af( i, n-k-l+j )
273 10 CONTINUE
274 20 CONTINUE
275*
276 IF( m-k-l.LT.0 ) THEN
277 DO 40 i = m + 1, k + l
278 DO 30 j = i, k + l
279 r( i, j ) = bf( i-k, n-k-l+j )
280 30 CONTINUE
281 40 CONTINUE
282 END IF
283*
284* Compute A:= U'*A*Q - D1*R
285*
286 CALL zgemm( 'No transpose', 'No transpose', m, n, n, cone, a, lda,
287 $ q, ldq, czero, work, lda )
288*
289 CALL zgemm( 'Conjugate transpose', 'No transpose', m, n, m, cone,
290 $ u, ldu, work, lda, czero, a, lda )
291*
292 DO 60 i = 1, k
293 DO 50 j = i, k + l
294 a( i, n-k-l+j ) = a( i, n-k-l+j ) - r( i, j )
295 50 CONTINUE
296 60 CONTINUE
297*
298 DO 80 i = k + 1, min( k+l, m )
299 DO 70 j = i, k + l
300 a( i, n-k-l+j ) = a( i, n-k-l+j ) - alpha( i )*r( i, j )
301 70 CONTINUE
302 80 CONTINUE
303*
304* Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) .
305*
306 resid = zlange( '1', m, n, a, lda, rwork )
307 IF( anorm.GT.zero ) THEN
308 result( 1 ) = ( ( resid / dble( max( 1, m, n ) ) ) / anorm ) /
309 $ ulp
310 ELSE
311 result( 1 ) = zero
312 END IF
313*
314* Compute B := V'*B*Q - D2*R
315*
316 CALL zgemm( 'No transpose', 'No transpose', p, n, n, cone, b, ldb,
317 $ q, ldq, czero, work, ldb )
318*
319 CALL zgemm( 'Conjugate transpose', 'No transpose', p, n, p, cone,
320 $ v, ldv, work, ldb, czero, b, ldb )
321*
322 DO 100 i = 1, l
323 DO 90 j = i, l
324 b( i, n-l+j ) = b( i, n-l+j ) - beta( k+i )*r( k+i, k+j )
325 90 CONTINUE
326 100 CONTINUE
327*
328* Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) .
329*
330 resid = zlange( '1', p, n, b, ldb, rwork )
331 IF( bnorm.GT.zero ) THEN
332 result( 2 ) = ( ( resid / dble( max( 1, p, n ) ) ) / bnorm ) /
333 $ ulp
334 ELSE
335 result( 2 ) = zero
336 END IF
337*
338* Compute I - U'*U
339*
340 CALL zlaset( 'Full', m, m, czero, cone, work, ldq )
341 CALL zherk( 'Upper', 'Conjugate transpose', m, m, -one, u, ldu,
342 $ one, work, ldu )
343*
344* Compute norm( I - U'*U ) / ( M * ULP ) .
345*
346 resid = zlanhe( '1', 'Upper', m, work, ldu, rwork )
347 result( 3 ) = ( resid / dble( max( 1, m ) ) ) / ulp
348*
349* Compute I - V'*V
350*
351 CALL zlaset( 'Full', p, p, czero, cone, work, ldv )
352 CALL zherk( 'Upper', 'Conjugate transpose', p, p, -one, v, ldv,
353 $ one, work, ldv )
354*
355* Compute norm( I - V'*V ) / ( P * ULP ) .
356*
357 resid = zlanhe( '1', 'Upper', p, work, ldv, rwork )
358 result( 4 ) = ( resid / dble( max( 1, p ) ) ) / ulp
359*
360* Compute I - Q'*Q
361*
362 CALL zlaset( 'Full', n, n, czero, cone, work, ldq )
363 CALL zherk( 'Upper', 'Conjugate transpose', n, n, -one, q, ldq,
364 $ one, work, ldq )
365*
366* Compute norm( I - Q'*Q ) / ( N * ULP ) .
367*
368 resid = zlanhe( '1', 'Upper', n, work, ldq, rwork )
369 result( 5 ) = ( resid / dble( max( 1, n ) ) ) / ulp
370*
371* Check sorting
372*
373 CALL dcopy( n, alpha, 1, rwork, 1 )
374 DO 110 i = k + 1, min( k+l, m )
375 j = iwork( i )
376 IF( i.NE.j ) THEN
377 temp = rwork( i )
378 rwork( i ) = rwork( j )
379 rwork( j ) = temp
380 END IF
381 110 CONTINUE
382*
383 result( 6 ) = zero
384 DO 120 i = k + 1, min( k+l, m ) - 1
385 IF( rwork( i ).LT.rwork( i+1 ) )
386 $ result( 6 ) = ulpinv
387 120 CONTINUE
388*
389 RETURN
390*
391* End of ZGSVTS3
392*
393 END
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
zggsvd3
subroutine zggsvd3(jobu, jobv, jobq, m, n, p, k, l, a, lda, b, ldb, alpha, beta, u, ldu, v, ldv, q, ldq, work, lwork, rwork, iwork, info)
ZGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices
Definition zggsvd3.f:351
zherk
subroutine zherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
ZHERK
Definition zherk.f:173
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:101
zlaset
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:104
zgsvts3
subroutine zgsvts3(m, p, n, a, af, lda, b, bf, ldb, u, ldu, v, ldv, q, ldq, alpha, beta, r, ldr, iwork, work, lwork, rwork, result)
ZGSVTS3
Definition zgsvts3.f:209