1 SUBROUTINE SGRQTS( M, P, N, A, AF, Q, R, LDA, TAUA, B, BF, Z, T,
2 $ BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER LDA, LDB, LWORK, M, P, N
10 * ..
11 * .. Array Arguments ..
12 REAL A( LDA, * ), AF( LDA, * ), R( LDA, * ),
13 $ Q( LDA, * ),
14 $ B( LDB, * ), BF( LDB, * ), T( LDB, * ),
15 $ Z( LDB, * ), BWK( LDB, * ),
16 $ TAUA( * ), TAUB( * ),
17 $ RESULT( 4 ), RWORK( * ), WORK( LWORK )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * SGRQTS tests SGGRQF, which computes the GRQ factorization of an
24 * M-by-N matrix A and a P-by-N matrix B: A = R*Q and B = Z*T*Q.
25 *
26 * Arguments
27 * =========
28 *
29 * M (input) INTEGER
30 * The number of rows of the matrix A. M >= 0.
31 *
32 * P (input) INTEGER
33 * The number of rows of the matrix B. P >= 0.
34 *
35 * N (input) INTEGER
36 * The number of columns of the matrices A and B. N >= 0.
37 *
38 * A (input) REAL array, dimension (LDA,N)
39 * The M-by-N matrix A.
40 *
41 * AF (output) REAL array, dimension (LDA,N)
42 * Details of the GRQ factorization of A and B, as returned
43 * by SGGRQF, see SGGRQF for further details.
44 *
45 * Q (output) REAL array, dimension (LDA,N)
46 * The N-by-N orthogonal matrix Q.
47 *
48 * R (workspace) REAL array, dimension (LDA,MAX(M,N))
49 *
50 * LDA (input) INTEGER
51 * The leading dimension of the arrays A, AF, R and Q.
52 * LDA >= max(M,N).
53 *
54 * TAUA (output) REAL array, dimension (min(M,N))
55 * The scalar factors of the elementary reflectors, as returned
56 * by SGGQRC.
57 *
58 * B (input) REAL array, dimension (LDB,N)
59 * On entry, the P-by-N matrix A.
60 *
61 * BF (output) REAL array, dimension (LDB,N)
62 * Details of the GQR factorization of A and B, as returned
63 * by SGGRQF, see SGGRQF for further details.
64 *
65 * Z (output) REAL array, dimension (LDB,P)
66 * The P-by-P orthogonal matrix Z.
67 *
68 * T (workspace) REAL array, dimension (LDB,max(P,N))
69 *
70 * BWK (workspace) REAL array, dimension (LDB,N)
71 *
72 * LDB (input) INTEGER
73 * The leading dimension of the arrays B, BF, Z and T.
74 * LDB >= max(P,N).
75 *
76 * TAUB (output) REAL array, dimension (min(P,N))
77 * The scalar factors of the elementary reflectors, as returned
78 * by SGGRQF.
79 *
80 * WORK (workspace) REAL array, dimension (LWORK)
81 *
82 * LWORK (input) INTEGER
83 * The dimension of the array WORK, LWORK >= max(M,P,N)**2.
84 *
85 * RWORK (workspace) REAL array, dimension (M)
86 *
87 * RESULT (output) REAL array, dimension (4)
88 * The test ratios:
89 * RESULT(1) = norm( R - A*Q' ) / ( MAX(M,N)*norm(A)*ULP)
90 * RESULT(2) = norm( T*Q - Z'*B ) / (MAX(P,N)*norm(B)*ULP)
91 * RESULT(3) = norm( I - Q'*Q ) / ( N*ULP )
92 * RESULT(4) = norm( I - Z'*Z ) / ( P*ULP )
93 *
94 * =====================================================================
95 *
96 * .. Parameters ..
97 REAL ZERO, ONE
98 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
99 REAL ROGUE
100 PARAMETER ( ROGUE = -1.0E+10 )
101 * ..
102 * .. Local Scalars ..
103 INTEGER INFO
104 REAL ANORM, BNORM, ULP, UNFL, RESID
105 * ..
106 * .. External Functions ..
107 REAL SLAMCH, SLANGE, SLANSY
108 EXTERNAL SLAMCH, SLANGE, SLANSY
109 * ..
110 * .. External Subroutines ..
111 EXTERNAL SGEMM, SGGRQF, SLACPY, SLASET, SORGQR,
112 $ SORGRQ, SSYRK
113 * ..
114 * .. Intrinsic Functions ..
115 INTRINSIC MAX, MIN, REAL
116 * ..
117 * .. Executable Statements ..
118 *
119 ULP = SLAMCH( 'Precision' )
120 UNFL = SLAMCH( 'Safe minimum' )
121 *
122 * Copy the matrix A to the array AF.
123 *
124 CALL SLACPY( 'Full', M, N, A, LDA, AF, LDA )
125 CALL SLACPY( 'Full', P, N, B, LDB, BF, LDB )
126 *
127 ANORM = MAX( SLANGE( '1', M, N, A, LDA, RWORK ), UNFL )
128 BNORM = MAX( SLANGE( '1', P, N, B, LDB, RWORK ), UNFL )
129 *
130 * Factorize the matrices A and B in the arrays AF and BF.
131 *
132 CALL SGGRQF( M, P, N, AF, LDA, TAUA, BF, LDB, TAUB, WORK,
133 $ LWORK, INFO )
134 *
135 * Generate the N-by-N matrix Q
136 *
137 CALL SLASET( 'Full', N, N, ROGUE, ROGUE, Q, LDA )
138 IF( M.LE.N ) THEN
139 IF( M.GT.0 .AND. M.LT.N )
140 $ CALL SLACPY( 'Full', M, N-M, AF, LDA, Q( N-M+1, 1 ), LDA )
141 IF( M.GT.1 )
142 $ CALL SLACPY( 'Lower', M-1, M-1, AF( 2, N-M+1 ), LDA,
143 $ Q( N-M+2, N-M+1 ), LDA )
144 ELSE
145 IF( N.GT.1 )
146 $ CALL SLACPY( 'Lower', N-1, N-1, AF( M-N+2, 1 ), LDA,
147 $ Q( 2, 1 ), LDA )
148 END IF
149 CALL SORGRQ( N, N, MIN( M, N ), Q, LDA, TAUA, WORK, LWORK, INFO )
150 *
151 * Generate the P-by-P matrix Z
152 *
153 CALL SLASET( 'Full', P, P, ROGUE, ROGUE, Z, LDB )
154 IF( P.GT.1 )
155 $ CALL SLACPY( 'Lower', P-1, N, BF( 2,1 ), LDB, Z( 2,1 ), LDB )
156 CALL SORGQR( P, P, MIN( P,N ), Z, LDB, TAUB, WORK, LWORK, INFO )
157 *
158 * Copy R
159 *
160 CALL SLASET( 'Full', M, N, ZERO, ZERO, R, LDA )
161 IF( M.LE.N )THEN
162 CALL SLACPY( 'Upper', M, M, AF( 1, N-M+1 ), LDA, R( 1, N-M+1 ),
163 $ LDA )
164 ELSE
165 CALL SLACPY( 'Full', M-N, N, AF, LDA, R, LDA )
166 CALL SLACPY( 'Upper', N, N, AF( M-N+1, 1 ), LDA, R( M-N+1, 1 ),
167 $ LDA )
168 END IF
169 *
170 * Copy T
171 *
172 CALL SLASET( 'Full', P, N, ZERO, ZERO, T, LDB )
173 CALL SLACPY( 'Upper', P, N, BF, LDB, T, LDB )
174 *
175 * Compute R - A*Q'
176 *
177 CALL SGEMM( 'No transpose', 'Transpose', M, N, N, -ONE, A, LDA, Q,
178 $ LDA, ONE, R, LDA )
179 *
180 * Compute norm( R - A*Q' ) / ( MAX(M,N)*norm(A)*ULP ) .
181 *
182 RESID = SLANGE( '1', M, N, R, LDA, RWORK )
183 IF( ANORM.GT.ZERO ) THEN
184 RESULT( 1 ) = ( ( RESID / REAL(MAX(1,M,N) ) ) / ANORM ) / ULP
185 ELSE
186 RESULT( 1 ) = ZERO
187 END IF
188 *
189 * Compute T*Q - Z'*B
190 *
191 CALL SGEMM( 'Transpose', 'No transpose', P, N, P, ONE, Z, LDB, B,
192 $ LDB, ZERO, BWK, LDB )
193 CALL SGEMM( 'No transpose', 'No transpose', P, N, N, ONE, T, LDB,
194 $ Q, LDA, -ONE, BWK, LDB )
195 *
196 * Compute norm( T*Q - Z'*B ) / ( MAX(P,N)*norm(A)*ULP ) .
197 *
198 RESID = SLANGE( '1', P, N, BWK, LDB, RWORK )
199 IF( BNORM.GT.ZERO ) THEN
200 RESULT( 2 ) = ( ( RESID / REAL( MAX( 1,P,M ) ) )/BNORM ) / ULP
201 ELSE
202 RESULT( 2 ) = ZERO
203 END IF
204 *
205 * Compute I - Q*Q'
206 *
207 CALL SLASET( 'Full', N, N, ZERO, ONE, R, LDA )
208 CALL SSYRK( 'Upper', 'No Transpose', N, N, -ONE, Q, LDA, ONE, R,
209 $ LDA )
210 *
211 * Compute norm( I - Q'*Q ) / ( N * ULP ) .
212 *
213 RESID = SLANSY( '1', 'Upper', N, R, LDA, RWORK )
214 RESULT( 3 ) = ( RESID / REAL( MAX( 1,N ) ) ) / ULP
215 *
216 * Compute I - Z'*Z
217 *
218 CALL SLASET( 'Full', P, P, ZERO, ONE, T, LDB )
219 CALL SSYRK( 'Upper', 'Transpose', P, P, -ONE, Z, LDB, ONE, T,
220 $ LDB )
221 *
222 * Compute norm( I - Z'*Z ) / ( P*ULP ) .
223 *
224 RESID = SLANSY( '1', 'Upper', P, T, LDB, RWORK )
225 RESULT( 4 ) = ( RESID / REAL( MAX( 1,P ) ) ) / ULP
226 *
227 RETURN
228 *
229 * End of SGRQTS
230 *
231 END
2 $ BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER LDA, LDB, LWORK, M, P, N
10 * ..
11 * .. Array Arguments ..
12 REAL A( LDA, * ), AF( LDA, * ), R( LDA, * ),
13 $ Q( LDA, * ),
14 $ B( LDB, * ), BF( LDB, * ), T( LDB, * ),
15 $ Z( LDB, * ), BWK( LDB, * ),
16 $ TAUA( * ), TAUB( * ),
17 $ RESULT( 4 ), RWORK( * ), WORK( LWORK )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * SGRQTS tests SGGRQF, which computes the GRQ factorization of an
24 * M-by-N matrix A and a P-by-N matrix B: A = R*Q and B = Z*T*Q.
25 *
26 * Arguments
27 * =========
28 *
29 * M (input) INTEGER
30 * The number of rows of the matrix A. M >= 0.
31 *
32 * P (input) INTEGER
33 * The number of rows of the matrix B. P >= 0.
34 *
35 * N (input) INTEGER
36 * The number of columns of the matrices A and B. N >= 0.
37 *
38 * A (input) REAL array, dimension (LDA,N)
39 * The M-by-N matrix A.
40 *
41 * AF (output) REAL array, dimension (LDA,N)
42 * Details of the GRQ factorization of A and B, as returned
43 * by SGGRQF, see SGGRQF for further details.
44 *
45 * Q (output) REAL array, dimension (LDA,N)
46 * The N-by-N orthogonal matrix Q.
47 *
48 * R (workspace) REAL array, dimension (LDA,MAX(M,N))
49 *
50 * LDA (input) INTEGER
51 * The leading dimension of the arrays A, AF, R and Q.
52 * LDA >= max(M,N).
53 *
54 * TAUA (output) REAL array, dimension (min(M,N))
55 * The scalar factors of the elementary reflectors, as returned
56 * by SGGQRC.
57 *
58 * B (input) REAL array, dimension (LDB,N)
59 * On entry, the P-by-N matrix A.
60 *
61 * BF (output) REAL array, dimension (LDB,N)
62 * Details of the GQR factorization of A and B, as returned
63 * by SGGRQF, see SGGRQF for further details.
64 *
65 * Z (output) REAL array, dimension (LDB,P)
66 * The P-by-P orthogonal matrix Z.
67 *
68 * T (workspace) REAL array, dimension (LDB,max(P,N))
69 *
70 * BWK (workspace) REAL array, dimension (LDB,N)
71 *
72 * LDB (input) INTEGER
73 * The leading dimension of the arrays B, BF, Z and T.
74 * LDB >= max(P,N).
75 *
76 * TAUB (output) REAL array, dimension (min(P,N))
77 * The scalar factors of the elementary reflectors, as returned
78 * by SGGRQF.
79 *
80 * WORK (workspace) REAL array, dimension (LWORK)
81 *
82 * LWORK (input) INTEGER
83 * The dimension of the array WORK, LWORK >= max(M,P,N)**2.
84 *
85 * RWORK (workspace) REAL array, dimension (M)
86 *
87 * RESULT (output) REAL array, dimension (4)
88 * The test ratios:
89 * RESULT(1) = norm( R - A*Q' ) / ( MAX(M,N)*norm(A)*ULP)
90 * RESULT(2) = norm( T*Q - Z'*B ) / (MAX(P,N)*norm(B)*ULP)
91 * RESULT(3) = norm( I - Q'*Q ) / ( N*ULP )
92 * RESULT(4) = norm( I - Z'*Z ) / ( P*ULP )
93 *
94 * =====================================================================
95 *
96 * .. Parameters ..
97 REAL ZERO, ONE
98 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
99 REAL ROGUE
100 PARAMETER ( ROGUE = -1.0E+10 )
101 * ..
102 * .. Local Scalars ..
103 INTEGER INFO
104 REAL ANORM, BNORM, ULP, UNFL, RESID
105 * ..
106 * .. External Functions ..
107 REAL SLAMCH, SLANGE, SLANSY
108 EXTERNAL SLAMCH, SLANGE, SLANSY
109 * ..
110 * .. External Subroutines ..
111 EXTERNAL SGEMM, SGGRQF, SLACPY, SLASET, SORGQR,
112 $ SORGRQ, SSYRK
113 * ..
114 * .. Intrinsic Functions ..
115 INTRINSIC MAX, MIN, REAL
116 * ..
117 * .. Executable Statements ..
118 *
119 ULP = SLAMCH( 'Precision' )
120 UNFL = SLAMCH( 'Safe minimum' )
121 *
122 * Copy the matrix A to the array AF.
123 *
124 CALL SLACPY( 'Full', M, N, A, LDA, AF, LDA )
125 CALL SLACPY( 'Full', P, N, B, LDB, BF, LDB )
126 *
127 ANORM = MAX( SLANGE( '1', M, N, A, LDA, RWORK ), UNFL )
128 BNORM = MAX( SLANGE( '1', P, N, B, LDB, RWORK ), UNFL )
129 *
130 * Factorize the matrices A and B in the arrays AF and BF.
131 *
132 CALL SGGRQF( M, P, N, AF, LDA, TAUA, BF, LDB, TAUB, WORK,
133 $ LWORK, INFO )
134 *
135 * Generate the N-by-N matrix Q
136 *
137 CALL SLASET( 'Full', N, N, ROGUE, ROGUE, Q, LDA )
138 IF( M.LE.N ) THEN
139 IF( M.GT.0 .AND. M.LT.N )
140 $ CALL SLACPY( 'Full', M, N-M, AF, LDA, Q( N-M+1, 1 ), LDA )
141 IF( M.GT.1 )
142 $ CALL SLACPY( 'Lower', M-1, M-1, AF( 2, N-M+1 ), LDA,
143 $ Q( N-M+2, N-M+1 ), LDA )
144 ELSE
145 IF( N.GT.1 )
146 $ CALL SLACPY( 'Lower', N-1, N-1, AF( M-N+2, 1 ), LDA,
147 $ Q( 2, 1 ), LDA )
148 END IF
149 CALL SORGRQ( N, N, MIN( M, N ), Q, LDA, TAUA, WORK, LWORK, INFO )
150 *
151 * Generate the P-by-P matrix Z
152 *
153 CALL SLASET( 'Full', P, P, ROGUE, ROGUE, Z, LDB )
154 IF( P.GT.1 )
155 $ CALL SLACPY( 'Lower', P-1, N, BF( 2,1 ), LDB, Z( 2,1 ), LDB )
156 CALL SORGQR( P, P, MIN( P,N ), Z, LDB, TAUB, WORK, LWORK, INFO )
157 *
158 * Copy R
159 *
160 CALL SLASET( 'Full', M, N, ZERO, ZERO, R, LDA )
161 IF( M.LE.N )THEN
162 CALL SLACPY( 'Upper', M, M, AF( 1, N-M+1 ), LDA, R( 1, N-M+1 ),
163 $ LDA )
164 ELSE
165 CALL SLACPY( 'Full', M-N, N, AF, LDA, R, LDA )
166 CALL SLACPY( 'Upper', N, N, AF( M-N+1, 1 ), LDA, R( M-N+1, 1 ),
167 $ LDA )
168 END IF
169 *
170 * Copy T
171 *
172 CALL SLASET( 'Full', P, N, ZERO, ZERO, T, LDB )
173 CALL SLACPY( 'Upper', P, N, BF, LDB, T, LDB )
174 *
175 * Compute R - A*Q'
176 *
177 CALL SGEMM( 'No transpose', 'Transpose', M, N, N, -ONE, A, LDA, Q,
178 $ LDA, ONE, R, LDA )
179 *
180 * Compute norm( R - A*Q' ) / ( MAX(M,N)*norm(A)*ULP ) .
181 *
182 RESID = SLANGE( '1', M, N, R, LDA, RWORK )
183 IF( ANORM.GT.ZERO ) THEN
184 RESULT( 1 ) = ( ( RESID / REAL(MAX(1,M,N) ) ) / ANORM ) / ULP
185 ELSE
186 RESULT( 1 ) = ZERO
187 END IF
188 *
189 * Compute T*Q - Z'*B
190 *
191 CALL SGEMM( 'Transpose', 'No transpose', P, N, P, ONE, Z, LDB, B,
192 $ LDB, ZERO, BWK, LDB )
193 CALL SGEMM( 'No transpose', 'No transpose', P, N, N, ONE, T, LDB,
194 $ Q, LDA, -ONE, BWK, LDB )
195 *
196 * Compute norm( T*Q - Z'*B ) / ( MAX(P,N)*norm(A)*ULP ) .
197 *
198 RESID = SLANGE( '1', P, N, BWK, LDB, RWORK )
199 IF( BNORM.GT.ZERO ) THEN
200 RESULT( 2 ) = ( ( RESID / REAL( MAX( 1,P,M ) ) )/BNORM ) / ULP
201 ELSE
202 RESULT( 2 ) = ZERO
203 END IF
204 *
205 * Compute I - Q*Q'
206 *
207 CALL SLASET( 'Full', N, N, ZERO, ONE, R, LDA )
208 CALL SSYRK( 'Upper', 'No Transpose', N, N, -ONE, Q, LDA, ONE, R,
209 $ LDA )
210 *
211 * Compute norm( I - Q'*Q ) / ( N * ULP ) .
212 *
213 RESID = SLANSY( '1', 'Upper', N, R, LDA, RWORK )
214 RESULT( 3 ) = ( RESID / REAL( MAX( 1,N ) ) ) / ULP
215 *
216 * Compute I - Z'*Z
217 *
218 CALL SLASET( 'Full', P, P, ZERO, ONE, T, LDB )
219 CALL SSYRK( 'Upper', 'Transpose', P, P, -ONE, Z, LDB, ONE, T,
220 $ LDB )
221 *
222 * Compute norm( I - Z'*Z ) / ( P*ULP ) .
223 *
224 RESID = SLANSY( '1', 'Upper', P, T, LDB, RWORK )
225 RESULT( 4 ) = ( RESID / REAL( MAX( 1,P ) ) ) / ULP
226 *
227 RETURN
228 *
229 * End of SGRQTS
230 *
231 END