LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 *> \date August 2015
204 *
205 *> \ingroup complex16_eig
206 *
207 * =====================================================================
208  SUBROUTINE zgsvts3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
209  $ ldv, q, ldq, alpha, beta, r, ldr, iwork, work,
210  $ lwork, rwork, result )
211 *
212 * -- LAPACK test routine (version 3.6.0) --
213 * -- LAPACK is a software package provided by Univ. of Tennessee, --
214 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
215 * August 2015
216 *
217 * .. Scalar Arguments ..
218  INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
219 * ..
220 * .. Array Arguments ..
221  INTEGER IWORK( * )
222  DOUBLE PRECISION ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
223  COMPLEX*16 A( lda, * ), AF( lda, * ), B( ldb, * ),
224  $ bf( ldb, * ), q( ldq, * ), r( ldr, * ),
225  $ u( ldu, * ), v( ldv, * ), work( lwork )
226 * ..
227 *
228 * =====================================================================
229 *
230 * .. Parameters ..
231  DOUBLE PRECISION ZERO, ONE
232  parameter ( zero = 0.0d+0, one = 1.0d+0 )
233  COMPLEX*16 CZERO, CONE
234  parameter ( czero = ( 0.0d+0, 0.0d+0 ),
235  $ cone = ( 1.0d+0, 0.0d+0 ) )
236 * ..
237 * .. Local Scalars ..
238  INTEGER I, INFO, J, K, L
239  DOUBLE PRECISION ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
240 * ..
241 * .. External Functions ..
242  DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE
243  EXTERNAL dlamch, zlange, zlanhe
244 * ..
245 * .. External Subroutines ..
246  EXTERNAL dcopy, zgemm, zggsvd3, zherk, zlacpy, zlaset
247 * ..
248 * .. Intrinsic Functions ..
249  INTRINSIC dble, max, min
250 * ..
251 * .. Executable Statements ..
252 *
253  ulp = dlamch( 'Precision' )
254  ulpinv = one / ulp
255  unfl = dlamch( 'Safe minimum' )
256 *
257 * Copy the matrix A to the array AF.
258 *
259  CALL zlacpy( 'Full', m, n, a, lda, af, lda )
260  CALL zlacpy( 'Full', p, n, b, ldb, bf, ldb )
261 *
262  anorm = max( zlange( '1', m, n, a, lda, rwork ), unfl )
263  bnorm = max( zlange( '1', p, n, b, ldb, rwork ), unfl )
264 *
265 * Factorize the matrices A and B in the arrays AF and BF.
266 *
267  CALL zggsvd3( 'U', 'V', 'Q', m, n, p, k, l, af, lda, bf, ldb,
268  $ alpha, beta, u, ldu, v, ldv, q, ldq, work, lwork,
269  $ rwork, iwork, info )
270 *
271 * Copy R
272 *
273  DO 20 i = 1, min( k+l, m )
274  DO 10 j = i, k + l
275  r( i, j ) = af( i, n-k-l+j )
276  10 CONTINUE
277  20 CONTINUE
278 *
279  IF( m-k-l.LT.0 ) THEN
280  DO 40 i = m + 1, k + l
281  DO 30 j = i, k + l
282  r( i, j ) = bf( i-k, n-k-l+j )
283  30 CONTINUE
284  40 CONTINUE
285  END IF
286 *
287 * Compute A:= U'*A*Q - D1*R
288 *
289  CALL zgemm( 'No transpose', 'No transpose', m, n, n, cone, a, lda,
290  $ q, ldq, czero, work, lda )
291 *
292  CALL zgemm( 'Conjugate transpose', 'No transpose', m, n, m, cone,
293  $ u, ldu, work, lda, czero, a, lda )
294 *
295  DO 60 i = 1, k
296  DO 50 j = i, k + l
297  a( i, n-k-l+j ) = a( i, n-k-l+j ) - r( i, j )
298  50 CONTINUE
299  60 CONTINUE
300 *
301  DO 80 i = k + 1, min( k+l, m )
302  DO 70 j = i, k + l
303  a( i, n-k-l+j ) = a( i, n-k-l+j ) - alpha( i )*r( i, j )
304  70 CONTINUE
305  80 CONTINUE
306 *
307 * Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) .
308 *
309  resid = zlange( '1', m, n, a, lda, rwork )
310  IF( anorm.GT.zero ) THEN
311  result( 1 ) = ( ( resid / dble( max( 1, m, n ) ) ) / anorm ) /
312  $ ulp
313  ELSE
314  result( 1 ) = zero
315  END IF
316 *
317 * Compute B := V'*B*Q - D2*R
318 *
319  CALL zgemm( 'No transpose', 'No transpose', p, n, n, cone, b, ldb,
320  $ q, ldq, czero, work, ldb )
321 *
322  CALL zgemm( 'Conjugate transpose', 'No transpose', p, n, p, cone,
323  $ v, ldv, work, ldb, czero, b, ldb )
324 *
325  DO 100 i = 1, l
326  DO 90 j = i, l
327  b( i, n-l+j ) = b( i, n-l+j ) - beta( k+i )*r( k+i, k+j )
328  90 CONTINUE
329  100 CONTINUE
330 *
331 * Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) .
332 *
333  resid = zlange( '1', p, n, b, ldb, rwork )
334  IF( bnorm.GT.zero ) THEN
335  result( 2 ) = ( ( resid / dble( max( 1, p, n ) ) ) / bnorm ) /
336  $ ulp
337  ELSE
338  result( 2 ) = zero
339  END IF
340 *
341 * Compute I - U'*U
342 *
343  CALL zlaset( 'Full', m, m, czero, cone, work, ldq )
344  CALL zherk( 'Upper', 'Conjugate transpose', m, m, -one, u, ldu,
345  $ one, work, ldu )
346 *
347 * Compute norm( I - U'*U ) / ( M * ULP ) .
348 *
349  resid = zlanhe( '1', 'Upper', m, work, ldu, rwork )
350  result( 3 ) = ( resid / dble( max( 1, m ) ) ) / ulp
351 *
352 * Compute I - V'*V
353 *
354  CALL zlaset( 'Full', p, p, czero, cone, work, ldv )
355  CALL zherk( 'Upper', 'Conjugate transpose', p, p, -one, v, ldv,
356  $ one, work, ldv )
357 *
358 * Compute norm( I - V'*V ) / ( P * ULP ) .
359 *
360  resid = zlanhe( '1', 'Upper', p, work, ldv, rwork )
361  result( 4 ) = ( resid / dble( max( 1, p ) ) ) / ulp
362 *
363 * Compute I - Q'*Q
364 *
365  CALL zlaset( 'Full', n, n, czero, cone, work, ldq )
366  CALL zherk( 'Upper', 'Conjugate transpose', n, n, -one, q, ldq,
367  $ one, work, ldq )
368 *
369 * Compute norm( I - Q'*Q ) / ( N * ULP ) .
370 *
371  resid = zlanhe( '1', 'Upper', n, work, ldq, rwork )
372  result( 5 ) = ( resid / dble( max( 1, n ) ) ) / ulp
373 *
374 * Check sorting
375 *
376  CALL dcopy( n, alpha, 1, rwork, 1 )
377  DO 110 i = k + 1, min( k+l, m )
378  j = iwork( i )
379  IF( i.NE.j ) THEN
380  temp = rwork( i )
381  rwork( i ) = rwork( j )
382  rwork( j ) = temp
383  END IF
384  110 CONTINUE
385 *
386  result( 6 ) = zero
387  DO 120 i = k + 1, min( k+l, m ) - 1
388  IF( rwork( i ).LT.rwork( i+1 ) )
389  $ result( 6 ) = ulpinv
390  120 CONTINUE
391 *
392  RETURN
393 *
394 * End of ZGSVTS3
395 *
396  END
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:105
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
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:355
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
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:108
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:211
zherk
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
Definition: zherk.f:175