1 SUBROUTINE CGSVTS( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
2 $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK,
3 $ LWORK, RWORK, RESULT )
4 *
5 * -- LAPACK test routine (version 3.1) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
11 * ..
12 * .. Array Arguments ..
13 INTEGER IWORK( * )
14 REAL ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
15 COMPLEX A( LDA, * ), AF( LDA, * ), B( LDB, * ),
16 $ BF( LDB, * ), Q( LDQ, * ), R( LDR, * ),
17 $ U( LDU, * ), V( LDV, * ), WORK( LWORK )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * CGSVTS tests CGGSVD, which computes the GSVD of an M-by-N matrix A
24 * and a P-by-N matrix B:
25 * U'*A*Q = D1*R and V'*B*Q = D2*R.
26 *
27 * Arguments
28 * =========
29 *
30 * M (input) INTEGER
31 * The number of rows of the matrix A. M >= 0.
32 *
33 * P (input) INTEGER
34 * The number of rows of the matrix B. P >= 0.
35 *
36 * N (input) INTEGER
37 * The number of columns of the matrices A and B. N >= 0.
38 *
39 * A (input) COMPLEX array, dimension (LDA,M)
40 * The M-by-N matrix A.
41 *
42 * AF (output) COMPLEX array, dimension (LDA,N)
43 * Details of the GSVD of A and B, as returned by CGGSVD,
44 * see CGGSVD for further details.
45 *
46 * LDA (input) INTEGER
47 * The leading dimension of the arrays A and AF.
48 * LDA >= max( 1,M ).
49 *
50 * B (input) COMPLEX array, dimension (LDB,P)
51 * On entry, the P-by-N matrix B.
52 *
53 * BF (output) COMPLEX array, dimension (LDB,N)
54 * Details of the GSVD of A and B, as returned by CGGSVD,
55 * see CGGSVD for further details.
56 *
57 * LDB (input) INTEGER
58 * The leading dimension of the arrays B and BF.
59 * LDB >= max(1,P).
60 *
61 * U (output) COMPLEX array, dimension(LDU,M)
62 * The M by M unitary matrix U.
63 *
64 * LDU (input) INTEGER
65 * The leading dimension of the array U. LDU >= max(1,M).
66 *
67 * V (output) COMPLEX array, dimension(LDV,M)
68 * The P by P unitary matrix V.
69 *
70 * LDV (input) INTEGER
71 * The leading dimension of the array V. LDV >= max(1,P).
72 *
73 * Q (output) COMPLEX array, dimension(LDQ,N)
74 * The N by N unitary matrix Q.
75 *
76 * LDQ (input) INTEGER
77 * The leading dimension of the array Q. LDQ >= max(1,N).
78 *
79 * ALPHA (output) REAL array, dimension (N)
80 * BETA (output) REAL array, dimension (N)
81 * The generalized singular value pairs of A and B, the
82 * ``diagonal'' matrices D1 and D2 are constructed from
83 * ALPHA and BETA, see subroutine CGGSVD for details.
84 *
85 * R (output) COMPLEX array, dimension(LDQ,N)
86 * The upper triangular matrix R.
87 *
88 * LDR (input) INTEGER
89 * The leading dimension of the array R. LDR >= max(1,N).
90 *
91 * IWORK (workspace) INTEGER array, dimension (N)
92 *
93 * WORK (workspace) COMPLEX array, dimension (LWORK)
94 *
95 * LWORK (input) INTEGER
96 * The dimension of the array WORK,
97 * LWORK >= max(M,P,N)*max(M,P,N).
98 *
99 * RWORK (workspace) REAL array, dimension (max(M,P,N))
100 *
101 * RESULT (output) REAL array, dimension (5)
102 * The test ratios:
103 * RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP)
104 * RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP)
105 * RESULT(3) = norm( I - U'*U ) / ( M*ULP )
106 * RESULT(4) = norm( I - V'*V ) / ( P*ULP )
107 * RESULT(5) = norm( I - Q'*Q ) / ( N*ULP )
108 * RESULT(6) = 0 if ALPHA is in decreasing order;
109 * = ULPINV otherwise.
110 *
111 * =====================================================================
112 *
113 * .. Parameters ..
114 REAL ZERO, ONE
115 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
116 COMPLEX CZERO, CONE
117 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
118 $ CONE = ( 1.0E+0, 0.0E+0 ) )
119 * ..
120 * .. Local Scalars ..
121 INTEGER I, INFO, J, K, L
122 REAL ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
123 * ..
124 * .. External Functions ..
125 REAL CLANGE, CLANHE, SLAMCH
126 EXTERNAL CLANGE, CLANHE, SLAMCH
127 * ..
128 * .. External Subroutines ..
129 EXTERNAL CGEMM, CGGSVD, CHERK, CLACPY, CLASET, SCOPY
130 * ..
131 * .. Intrinsic Functions ..
132 INTRINSIC MAX, MIN, REAL
133 * ..
134 * .. Executable Statements ..
135 *
136 ULP = SLAMCH( 'Precision' )
137 ULPINV = ONE / ULP
138 UNFL = SLAMCH( 'Safe minimum' )
139 *
140 * Copy the matrix A to the array AF.
141 *
142 CALL CLACPY( 'Full', M, N, A, LDA, AF, LDA )
143 CALL CLACPY( 'Full', P, N, B, LDB, BF, LDB )
144 *
145 ANORM = MAX( CLANGE( '1', M, N, A, LDA, RWORK ), UNFL )
146 BNORM = MAX( CLANGE( '1', P, N, B, LDB, RWORK ), UNFL )
147 *
148 * Factorize the matrices A and B in the arrays AF and BF.
149 *
150 CALL CGGSVD( 'U', 'V', 'Q', M, N, P, K, L, AF, LDA, BF, LDB,
151 $ ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, RWORK,
152 $ IWORK, INFO )
153 *
154 * Copy R
155 *
156 DO 20 I = 1, MIN( K+L, M )
157 DO 10 J = I, K + L
158 R( I, J ) = AF( I, N-K-L+J )
159 10 CONTINUE
160 20 CONTINUE
161 *
162 IF( M-K-L.LT.0 ) THEN
163 DO 40 I = M + 1, K + L
164 DO 30 J = I, K + L
165 R( I, J ) = BF( I-K, N-K-L+J )
166 30 CONTINUE
167 40 CONTINUE
168 END IF
169 *
170 * Compute A:= U'*A*Q - D1*R
171 *
172 CALL CGEMM( 'No transpose', 'No transpose', M, N, N, CONE, A, LDA,
173 $ Q, LDQ, CZERO, WORK, LDA )
174 *
175 CALL CGEMM( 'Conjugate transpose', 'No transpose', M, N, M, CONE,
176 $ U, LDU, WORK, LDA, CZERO, A, LDA )
177 *
178 DO 60 I = 1, K
179 DO 50 J = I, K + L
180 A( I, N-K-L+J ) = A( I, N-K-L+J ) - R( I, J )
181 50 CONTINUE
182 60 CONTINUE
183 *
184 DO 80 I = K + 1, MIN( K+L, M )
185 DO 70 J = I, K + L
186 A( I, N-K-L+J ) = A( I, N-K-L+J ) - ALPHA( I )*R( I, J )
187 70 CONTINUE
188 80 CONTINUE
189 *
190 * Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) .
191 *
192 RESID = CLANGE( '1', M, N, A, LDA, RWORK )
193 IF( ANORM.GT.ZERO ) THEN
194 RESULT( 1 ) = ( ( RESID / REAL( MAX( 1, M, N ) ) ) / ANORM ) /
195 $ ULP
196 ELSE
197 RESULT( 1 ) = ZERO
198 END IF
199 *
200 * Compute B := V'*B*Q - D2*R
201 *
202 CALL CGEMM( 'No transpose', 'No transpose', P, N, N, CONE, B, LDB,
203 $ Q, LDQ, CZERO, WORK, LDB )
204 *
205 CALL CGEMM( 'Conjugate transpose', 'No transpose', P, N, P, CONE,
206 $ V, LDV, WORK, LDB, CZERO, B, LDB )
207 *
208 DO 100 I = 1, L
209 DO 90 J = I, L
210 B( I, N-L+J ) = B( I, N-L+J ) - BETA( K+I )*R( K+I, K+J )
211 90 CONTINUE
212 100 CONTINUE
213 *
214 * Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) .
215 *
216 RESID = CLANGE( '1', P, N, B, LDB, RWORK )
217 IF( BNORM.GT.ZERO ) THEN
218 RESULT( 2 ) = ( ( RESID / REAL( MAX( 1, P, N ) ) ) / BNORM ) /
219 $ ULP
220 ELSE
221 RESULT( 2 ) = ZERO
222 END IF
223 *
224 * Compute I - U'*U
225 *
226 CALL CLASET( 'Full', M, M, CZERO, CONE, WORK, LDQ )
227 CALL CHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, U, LDU,
228 $ ONE, WORK, LDU )
229 *
230 * Compute norm( I - U'*U ) / ( M * ULP ) .
231 *
232 RESID = CLANHE( '1', 'Upper', M, WORK, LDU, RWORK )
233 RESULT( 3 ) = ( RESID / REAL( MAX( 1, M ) ) ) / ULP
234 *
235 * Compute I - V'*V
236 *
237 CALL CLASET( 'Full', P, P, CZERO, CONE, WORK, LDV )
238 CALL CHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, V, LDV,
239 $ ONE, WORK, LDV )
240 *
241 * Compute norm( I - V'*V ) / ( P * ULP ) .
242 *
243 RESID = CLANHE( '1', 'Upper', P, WORK, LDV, RWORK )
244 RESULT( 4 ) = ( RESID / REAL( MAX( 1, P ) ) ) / ULP
245 *
246 * Compute I - Q'*Q
247 *
248 CALL CLASET( 'Full', N, N, CZERO, CONE, WORK, LDQ )
249 CALL CHERK( 'Upper', 'Conjugate transpose', N, N, -ONE, Q, LDQ,
250 $ ONE, WORK, LDQ )
251 *
252 * Compute norm( I - Q'*Q ) / ( N * ULP ) .
253 *
254 RESID = CLANHE( '1', 'Upper', N, WORK, LDQ, RWORK )
255 RESULT( 5 ) = ( RESID / REAL( MAX( 1, N ) ) ) / ULP
256 *
257 * Check sorting
258 *
259 CALL SCOPY( N, ALPHA, 1, RWORK, 1 )
260 DO 110 I = K + 1, MIN( K+L, M )
261 J = IWORK( I )
262 IF( I.NE.J ) THEN
263 TEMP = RWORK( I )
264 RWORK( I ) = RWORK( J )
265 RWORK( J ) = TEMP
266 END IF
267 110 CONTINUE
268 *
269 RESULT( 6 ) = ZERO
270 DO 120 I = K + 1, MIN( K+L, M ) - 1
271 IF( RWORK( I ).LT.RWORK( I+1 ) )
272 $ RESULT( 6 ) = ULPINV
273 120 CONTINUE
274 *
275 RETURN
276 *
277 * End of CGSVTS
278 *
279 END
2 $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK,
3 $ LWORK, RWORK, RESULT )
4 *
5 * -- LAPACK test routine (version 3.1) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
11 * ..
12 * .. Array Arguments ..
13 INTEGER IWORK( * )
14 REAL ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * )
15 COMPLEX A( LDA, * ), AF( LDA, * ), B( LDB, * ),
16 $ BF( LDB, * ), Q( LDQ, * ), R( LDR, * ),
17 $ U( LDU, * ), V( LDV, * ), WORK( LWORK )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * CGSVTS tests CGGSVD, which computes the GSVD of an M-by-N matrix A
24 * and a P-by-N matrix B:
25 * U'*A*Q = D1*R and V'*B*Q = D2*R.
26 *
27 * Arguments
28 * =========
29 *
30 * M (input) INTEGER
31 * The number of rows of the matrix A. M >= 0.
32 *
33 * P (input) INTEGER
34 * The number of rows of the matrix B. P >= 0.
35 *
36 * N (input) INTEGER
37 * The number of columns of the matrices A and B. N >= 0.
38 *
39 * A (input) COMPLEX array, dimension (LDA,M)
40 * The M-by-N matrix A.
41 *
42 * AF (output) COMPLEX array, dimension (LDA,N)
43 * Details of the GSVD of A and B, as returned by CGGSVD,
44 * see CGGSVD for further details.
45 *
46 * LDA (input) INTEGER
47 * The leading dimension of the arrays A and AF.
48 * LDA >= max( 1,M ).
49 *
50 * B (input) COMPLEX array, dimension (LDB,P)
51 * On entry, the P-by-N matrix B.
52 *
53 * BF (output) COMPLEX array, dimension (LDB,N)
54 * Details of the GSVD of A and B, as returned by CGGSVD,
55 * see CGGSVD for further details.
56 *
57 * LDB (input) INTEGER
58 * The leading dimension of the arrays B and BF.
59 * LDB >= max(1,P).
60 *
61 * U (output) COMPLEX array, dimension(LDU,M)
62 * The M by M unitary matrix U.
63 *
64 * LDU (input) INTEGER
65 * The leading dimension of the array U. LDU >= max(1,M).
66 *
67 * V (output) COMPLEX array, dimension(LDV,M)
68 * The P by P unitary matrix V.
69 *
70 * LDV (input) INTEGER
71 * The leading dimension of the array V. LDV >= max(1,P).
72 *
73 * Q (output) COMPLEX array, dimension(LDQ,N)
74 * The N by N unitary matrix Q.
75 *
76 * LDQ (input) INTEGER
77 * The leading dimension of the array Q. LDQ >= max(1,N).
78 *
79 * ALPHA (output) REAL array, dimension (N)
80 * BETA (output) REAL array, dimension (N)
81 * The generalized singular value pairs of A and B, the
82 * ``diagonal'' matrices D1 and D2 are constructed from
83 * ALPHA and BETA, see subroutine CGGSVD for details.
84 *
85 * R (output) COMPLEX array, dimension(LDQ,N)
86 * The upper triangular matrix R.
87 *
88 * LDR (input) INTEGER
89 * The leading dimension of the array R. LDR >= max(1,N).
90 *
91 * IWORK (workspace) INTEGER array, dimension (N)
92 *
93 * WORK (workspace) COMPLEX array, dimension (LWORK)
94 *
95 * LWORK (input) INTEGER
96 * The dimension of the array WORK,
97 * LWORK >= max(M,P,N)*max(M,P,N).
98 *
99 * RWORK (workspace) REAL array, dimension (max(M,P,N))
100 *
101 * RESULT (output) REAL array, dimension (5)
102 * The test ratios:
103 * RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP)
104 * RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP)
105 * RESULT(3) = norm( I - U'*U ) / ( M*ULP )
106 * RESULT(4) = norm( I - V'*V ) / ( P*ULP )
107 * RESULT(5) = norm( I - Q'*Q ) / ( N*ULP )
108 * RESULT(6) = 0 if ALPHA is in decreasing order;
109 * = ULPINV otherwise.
110 *
111 * =====================================================================
112 *
113 * .. Parameters ..
114 REAL ZERO, ONE
115 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
116 COMPLEX CZERO, CONE
117 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
118 $ CONE = ( 1.0E+0, 0.0E+0 ) )
119 * ..
120 * .. Local Scalars ..
121 INTEGER I, INFO, J, K, L
122 REAL ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL
123 * ..
124 * .. External Functions ..
125 REAL CLANGE, CLANHE, SLAMCH
126 EXTERNAL CLANGE, CLANHE, SLAMCH
127 * ..
128 * .. External Subroutines ..
129 EXTERNAL CGEMM, CGGSVD, CHERK, CLACPY, CLASET, SCOPY
130 * ..
131 * .. Intrinsic Functions ..
132 INTRINSIC MAX, MIN, REAL
133 * ..
134 * .. Executable Statements ..
135 *
136 ULP = SLAMCH( 'Precision' )
137 ULPINV = ONE / ULP
138 UNFL = SLAMCH( 'Safe minimum' )
139 *
140 * Copy the matrix A to the array AF.
141 *
142 CALL CLACPY( 'Full', M, N, A, LDA, AF, LDA )
143 CALL CLACPY( 'Full', P, N, B, LDB, BF, LDB )
144 *
145 ANORM = MAX( CLANGE( '1', M, N, A, LDA, RWORK ), UNFL )
146 BNORM = MAX( CLANGE( '1', P, N, B, LDB, RWORK ), UNFL )
147 *
148 * Factorize the matrices A and B in the arrays AF and BF.
149 *
150 CALL CGGSVD( 'U', 'V', 'Q', M, N, P, K, L, AF, LDA, BF, LDB,
151 $ ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, RWORK,
152 $ IWORK, INFO )
153 *
154 * Copy R
155 *
156 DO 20 I = 1, MIN( K+L, M )
157 DO 10 J = I, K + L
158 R( I, J ) = AF( I, N-K-L+J )
159 10 CONTINUE
160 20 CONTINUE
161 *
162 IF( M-K-L.LT.0 ) THEN
163 DO 40 I = M + 1, K + L
164 DO 30 J = I, K + L
165 R( I, J ) = BF( I-K, N-K-L+J )
166 30 CONTINUE
167 40 CONTINUE
168 END IF
169 *
170 * Compute A:= U'*A*Q - D1*R
171 *
172 CALL CGEMM( 'No transpose', 'No transpose', M, N, N, CONE, A, LDA,
173 $ Q, LDQ, CZERO, WORK, LDA )
174 *
175 CALL CGEMM( 'Conjugate transpose', 'No transpose', M, N, M, CONE,
176 $ U, LDU, WORK, LDA, CZERO, A, LDA )
177 *
178 DO 60 I = 1, K
179 DO 50 J = I, K + L
180 A( I, N-K-L+J ) = A( I, N-K-L+J ) - R( I, J )
181 50 CONTINUE
182 60 CONTINUE
183 *
184 DO 80 I = K + 1, MIN( K+L, M )
185 DO 70 J = I, K + L
186 A( I, N-K-L+J ) = A( I, N-K-L+J ) - ALPHA( I )*R( I, J )
187 70 CONTINUE
188 80 CONTINUE
189 *
190 * Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) .
191 *
192 RESID = CLANGE( '1', M, N, A, LDA, RWORK )
193 IF( ANORM.GT.ZERO ) THEN
194 RESULT( 1 ) = ( ( RESID / REAL( MAX( 1, M, N ) ) ) / ANORM ) /
195 $ ULP
196 ELSE
197 RESULT( 1 ) = ZERO
198 END IF
199 *
200 * Compute B := V'*B*Q - D2*R
201 *
202 CALL CGEMM( 'No transpose', 'No transpose', P, N, N, CONE, B, LDB,
203 $ Q, LDQ, CZERO, WORK, LDB )
204 *
205 CALL CGEMM( 'Conjugate transpose', 'No transpose', P, N, P, CONE,
206 $ V, LDV, WORK, LDB, CZERO, B, LDB )
207 *
208 DO 100 I = 1, L
209 DO 90 J = I, L
210 B( I, N-L+J ) = B( I, N-L+J ) - BETA( K+I )*R( K+I, K+J )
211 90 CONTINUE
212 100 CONTINUE
213 *
214 * Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) .
215 *
216 RESID = CLANGE( '1', P, N, B, LDB, RWORK )
217 IF( BNORM.GT.ZERO ) THEN
218 RESULT( 2 ) = ( ( RESID / REAL( MAX( 1, P, N ) ) ) / BNORM ) /
219 $ ULP
220 ELSE
221 RESULT( 2 ) = ZERO
222 END IF
223 *
224 * Compute I - U'*U
225 *
226 CALL CLASET( 'Full', M, M, CZERO, CONE, WORK, LDQ )
227 CALL CHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, U, LDU,
228 $ ONE, WORK, LDU )
229 *
230 * Compute norm( I - U'*U ) / ( M * ULP ) .
231 *
232 RESID = CLANHE( '1', 'Upper', M, WORK, LDU, RWORK )
233 RESULT( 3 ) = ( RESID / REAL( MAX( 1, M ) ) ) / ULP
234 *
235 * Compute I - V'*V
236 *
237 CALL CLASET( 'Full', P, P, CZERO, CONE, WORK, LDV )
238 CALL CHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, V, LDV,
239 $ ONE, WORK, LDV )
240 *
241 * Compute norm( I - V'*V ) / ( P * ULP ) .
242 *
243 RESID = CLANHE( '1', 'Upper', P, WORK, LDV, RWORK )
244 RESULT( 4 ) = ( RESID / REAL( MAX( 1, P ) ) ) / ULP
245 *
246 * Compute I - Q'*Q
247 *
248 CALL CLASET( 'Full', N, N, CZERO, CONE, WORK, LDQ )
249 CALL CHERK( 'Upper', 'Conjugate transpose', N, N, -ONE, Q, LDQ,
250 $ ONE, WORK, LDQ )
251 *
252 * Compute norm( I - Q'*Q ) / ( N * ULP ) .
253 *
254 RESID = CLANHE( '1', 'Upper', N, WORK, LDQ, RWORK )
255 RESULT( 5 ) = ( RESID / REAL( MAX( 1, N ) ) ) / ULP
256 *
257 * Check sorting
258 *
259 CALL SCOPY( N, ALPHA, 1, RWORK, 1 )
260 DO 110 I = K + 1, MIN( K+L, M )
261 J = IWORK( I )
262 IF( I.NE.J ) THEN
263 TEMP = RWORK( I )
264 RWORK( I ) = RWORK( J )
265 RWORK( J ) = TEMP
266 END IF
267 110 CONTINUE
268 *
269 RESULT( 6 ) = ZERO
270 DO 120 I = K + 1, MIN( K+L, M ) - 1
271 IF( RWORK( I ).LT.RWORK( I+1 ) )
272 $ RESULT( 6 ) = ULPINV
273 120 CONTINUE
274 *
275 RETURN
276 *
277 * End of CGSVTS
278 *
279 END