LAPACK  3.5.0 LAPACK: Linear Algebra PACKage
cgsvts.f
Go to the documentation of this file.
1 *> \brief \b CGSVTS
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 CGSVTS( 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 * REAL ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
21 * COMPLEX 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 *> CGSVTS tests CGGSVD, 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 array, dimension (LDA,M)
61 *> The M-by-N matrix A.
62 *> \endverbatim
63 *>
64 *> \param[out] AF
65 *> \verbatim
66 *> AF is COMPLEX array, dimension (LDA,N)
67 *> Details of the GSVD of A and B, as returned by CGGSVD,
68 *> see CGGSVD 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 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 array, dimension (LDB,N)
87 *> Details of the GSVD of A and B, as returned by CGGSVD,
88 *> see CGGSVD 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 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 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 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 REAL array, dimension (N)
137 *> \endverbatim
138 *>
139 *> \param[out] BETA
140 *> \verbatim
141 *> BETA is REAL 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 CGGSVD for details.
146 *> \endverbatim
147 *>
148 *> \param[out] R
149 *> \verbatim
150 *> R is COMPLEX 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 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 REAL array, dimension (max(M,P,N))
180 *> \endverbatim
181 *>
182 *> \param[out] RESULT
183 *> \verbatim
184 *> RESULT is REAL array, dimension (5)
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 November 2011
204 *
205 *> \ingroup complex_eig
206 *
207 * =====================================================================
208  SUBROUTINE cgsvts( 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.4.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 * November 2011
216 *
217 * .. Scalar Arguments ..
218  INTEGER lda, ldb, ldq, ldr, ldu, ldv, lwork, m, n, p
219 * ..
220 * .. Array Arguments ..
221  INTEGER iwork( * )
222  REAL alpha( * ), beta( * ), result( 6 ), rwork( * )
223  COMPLEX 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  REAL zero, one
232  parameter( zero = 0.0e+0, one = 1.0e+0 )
233  COMPLEX czero, cone
234  parameter( czero = ( 0.0e+0, 0.0e+0 ),
235  \$ cone = ( 1.0e+0, 0.0e+0 ) )
236 * ..
237 * .. Local Scalars ..
238  INTEGER i, info, j, k, l
239  REAL anorm, bnorm, resid, temp, ulp, ulpinv, unfl
240 * ..
241 * .. External Functions ..
242  REAL clange, clanhe, slamch
243  EXTERNAL clange, clanhe, slamch
244 * ..
245 * .. External Subroutines ..
246  EXTERNAL cgemm, cggsvd, cherk, clacpy, claset, scopy
247 * ..
248 * .. Intrinsic Functions ..
249  INTRINSIC max, min, real
250 * ..
251 * .. Executable Statements ..
252 *
253  ulp = slamch( 'Precision' )
254  ulpinv = one / ulp
255  unfl = slamch( 'Safe minimum' )
256 *
257 * Copy the matrix A to the array AF.
258 *
259  CALL clacpy( 'Full', m, n, a, lda, af, lda )
260  CALL clacpy( 'Full', p, n, b, ldb, bf, ldb )
261 *
262  anorm = max( clange( '1', m, n, a, lda, rwork ), unfl )
263  bnorm = max( clange( '1', p, n, b, ldb, rwork ), unfl )
264 *
265 * Factorize the matrices A and B in the arrays AF and BF.
266 *
267  CALL cggsvd( 'U', 'V', 'Q', m, n, p, k, l, af, lda, bf, ldb,
268  \$ alpha, beta, u, ldu, v, ldv, q, ldq, work, rwork,
269  \$ 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 cgemm( 'No transpose', 'No transpose', m, n, n, cone, a, lda,
290  \$ q, ldq, czero, work, lda )
291 *
292  CALL cgemm( '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 = clange( '1', m, n, a, lda, rwork )
310  IF( anorm.GT.zero ) THEN
311  result( 1 ) = ( ( resid / REAL( 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 cgemm( 'No transpose', 'No transpose', p, n, n, cone, b, ldb,
320  \$ q, ldq, czero, work, ldb )
321 *
322  CALL cgemm( '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 = clange( '1', p, n, b, ldb, rwork )
334  IF( bnorm.GT.zero ) THEN
335  result( 2 ) = ( ( resid / REAL( 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 claset( 'Full', m, m, czero, cone, work, ldq )
344  CALL cherk( 'Upper', 'Conjugate transpose', m, m, -one, u, ldu,
345  \$ one, work, ldu )
346 *
347 * Compute norm( I - U'*U ) / ( M * ULP ) .
348 *
349  resid = clanhe( '1', 'Upper', m, work, ldu, rwork )
350  result( 3 ) = ( resid / REAL( MAX( 1, M ) ) ) / ulp
351 *
352 * Compute I - V'*V
353 *
354  CALL claset( 'Full', p, p, czero, cone, work, ldv )
355  CALL cherk( 'Upper', 'Conjugate transpose', p, p, -one, v, ldv,
356  \$ one, work, ldv )
357 *
358 * Compute norm( I - V'*V ) / ( P * ULP ) .
359 *
360  resid = clanhe( '1', 'Upper', p, work, ldv, rwork )
361  result( 4 ) = ( resid / REAL( MAX( 1, P ) ) ) / ulp
362 *
363 * Compute I - Q'*Q
364 *
365  CALL claset( 'Full', n, n, czero, cone, work, ldq )
366  CALL cherk( 'Upper', 'Conjugate transpose', n, n, -one, q, ldq,
367  \$ one, work, ldq )
368 *
369 * Compute norm( I - Q'*Q ) / ( N * ULP ) .
370 *
371  resid = clanhe( '1', 'Upper', n, work, ldq, rwork )
372  result( 5 ) = ( resid / REAL( MAX( 1, N ) ) ) / ulp
373 *
374 * Check sorting
375 *
376  CALL scopy( 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 CGSVTS
395 *
396  END