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