1 SUBROUTINE CGET51( ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK,
2 $ 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 ITYPE, LDA, LDB, LDU, LDV, N
10 REAL RESULT
11 * ..
12 * .. Array Arguments ..
13 REAL RWORK( * )
14 COMPLEX A( LDA, * ), B( LDB, * ), U( LDU, * ),
15 $ V( LDV, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * CGET51 generally checks a decomposition of the form
22 *
23 * A = U B V*
24 *
25 * where * means conjugate transpose and U and V are unitary.
26 *
27 * Specifically, if ITYPE=1
28 *
29 * RESULT = | A - U B V* | / ( |A| n ulp )
30 *
31 * If ITYPE=2, then:
32 *
33 * RESULT = | A - B | / ( |A| n ulp )
34 *
35 * If ITYPE=3, then:
36 *
37 * RESULT = | I - UU* | / ( n ulp )
38 *
39 * Arguments
40 * =========
41 *
42 * ITYPE (input) INTEGER
43 * Specifies the type of tests to be performed.
44 * =1: RESULT = | A - U B V* | / ( |A| n ulp )
45 * =2: RESULT = | A - B | / ( |A| n ulp )
46 * =3: RESULT = | I - UU* | / ( n ulp )
47 *
48 * N (input) INTEGER
49 * The size of the matrix. If it is zero, CGET51 does nothing.
50 * It must be at least zero.
51 *
52 * A (input) COMPLEX array, dimension (LDA, N)
53 * The original (unfactored) matrix.
54 *
55 * LDA (input) INTEGER
56 * The leading dimension of A. It must be at least 1
57 * and at least N.
58 *
59 * B (input) COMPLEX array, dimension (LDB, N)
60 * The factored matrix.
61 *
62 * LDB (input) INTEGER
63 * The leading dimension of B. It must be at least 1
64 * and at least N.
65 *
66 * U (input) COMPLEX array, dimension (LDU, N)
67 * The unitary matrix on the left-hand side in the
68 * decomposition.
69 * Not referenced if ITYPE=2
70 *
71 * LDU (input) INTEGER
72 * The leading dimension of U. LDU must be at least N and
73 * at least 1.
74 *
75 * V (input) COMPLEX array, dimension (LDV, N)
76 * The unitary matrix on the left-hand side in the
77 * decomposition.
78 * Not referenced if ITYPE=2
79 *
80 * LDV (input) INTEGER
81 * The leading dimension of V. LDV must be at least N and
82 * at least 1.
83 *
84 * WORK (workspace) COMPLEX array, dimension (2*N**2)
85 *
86 * RWORK (workspace) REAL array, dimension (N)
87 *
88 * RESULT (output) REAL
89 * The values computed by the test specified by ITYPE. The
90 * value is currently limited to 1/ulp, to avoid overflow.
91 * Errors are flagged by RESULT=10/ulp.
92 *
93 * =====================================================================
94 *
95 * .. Parameters ..
96 REAL ZERO, ONE, TEN
97 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 10.0E+0 )
98 COMPLEX CZERO, CONE
99 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
100 $ CONE = ( 1.0E+0, 0.0E+0 ) )
101 * ..
102 * .. Local Scalars ..
103 INTEGER JCOL, JDIAG, JROW
104 REAL ANORM, ULP, UNFL, WNORM
105 * ..
106 * .. External Functions ..
107 REAL CLANGE, SLAMCH
108 EXTERNAL CLANGE, SLAMCH
109 * ..
110 * .. External Subroutines ..
111 EXTERNAL CGEMM, CLACPY
112 * ..
113 * .. Intrinsic Functions ..
114 INTRINSIC MAX, MIN, REAL
115 * ..
116 * .. Executable Statements ..
117 *
118 RESULT = ZERO
119 IF( N.LE.0 )
120 $ RETURN
121 *
122 * Constants
123 *
124 UNFL = SLAMCH( 'Safe minimum' )
125 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
126 *
127 * Some Error Checks
128 *
129 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
130 RESULT = TEN / ULP
131 RETURN
132 END IF
133 *
134 IF( ITYPE.LE.2 ) THEN
135 *
136 * Tests scaled by the norm(A)
137 *
138 ANORM = MAX( CLANGE( '1', N, N, A, LDA, RWORK ), UNFL )
139 *
140 IF( ITYPE.EQ.1 ) THEN
141 *
142 * ITYPE=1: Compute W = A - UBV'
143 *
144 CALL CLACPY( ' ', N, N, A, LDA, WORK, N )
145 CALL CGEMM( 'N', 'N', N, N, N, CONE, U, LDU, B, LDB, CZERO,
146 $ WORK( N**2+1 ), N )
147 *
148 CALL CGEMM( 'N', 'C', N, N, N, -CONE, WORK( N**2+1 ), N, V,
149 $ LDV, CONE, WORK, N )
150 *
151 ELSE
152 *
153 * ITYPE=2: Compute W = A - B
154 *
155 CALL CLACPY( ' ', N, N, B, LDB, WORK, N )
156 *
157 DO 20 JCOL = 1, N
158 DO 10 JROW = 1, N
159 WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) )
160 $ - A( JROW, JCOL )
161 10 CONTINUE
162 20 CONTINUE
163 END IF
164 *
165 * Compute norm(W)/ ( ulp*norm(A) )
166 *
167 WNORM = CLANGE( '1', N, N, WORK, N, RWORK )
168 *
169 IF( ANORM.GT.WNORM ) THEN
170 RESULT = ( WNORM / ANORM ) / ( N*ULP )
171 ELSE
172 IF( ANORM.LT.ONE ) THEN
173 RESULT = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
174 ELSE
175 RESULT = MIN( WNORM / ANORM, REAL( N ) ) / ( N*ULP )
176 END IF
177 END IF
178 *
179 ELSE
180 *
181 * Tests not scaled by norm(A)
182 *
183 * ITYPE=3: Compute UU' - I
184 *
185 CALL CGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO,
186 $ WORK, N )
187 *
188 DO 30 JDIAG = 1, N
189 WORK( ( N+1 )*( JDIAG-1 )+1 ) = WORK( ( N+1 )*( JDIAG-1 )+
190 $ 1 ) - CONE
191 30 CONTINUE
192 *
193 RESULT = MIN( CLANGE( '1', N, N, WORK, N, RWORK ),
194 $ REAL( N ) ) / ( N*ULP )
195 END IF
196 *
197 RETURN
198 *
199 * End of CGET51
200 *
201 END
2 $ 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 ITYPE, LDA, LDB, LDU, LDV, N
10 REAL RESULT
11 * ..
12 * .. Array Arguments ..
13 REAL RWORK( * )
14 COMPLEX A( LDA, * ), B( LDB, * ), U( LDU, * ),
15 $ V( LDV, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * CGET51 generally checks a decomposition of the form
22 *
23 * A = U B V*
24 *
25 * where * means conjugate transpose and U and V are unitary.
26 *
27 * Specifically, if ITYPE=1
28 *
29 * RESULT = | A - U B V* | / ( |A| n ulp )
30 *
31 * If ITYPE=2, then:
32 *
33 * RESULT = | A - B | / ( |A| n ulp )
34 *
35 * If ITYPE=3, then:
36 *
37 * RESULT = | I - UU* | / ( n ulp )
38 *
39 * Arguments
40 * =========
41 *
42 * ITYPE (input) INTEGER
43 * Specifies the type of tests to be performed.
44 * =1: RESULT = | A - U B V* | / ( |A| n ulp )
45 * =2: RESULT = | A - B | / ( |A| n ulp )
46 * =3: RESULT = | I - UU* | / ( n ulp )
47 *
48 * N (input) INTEGER
49 * The size of the matrix. If it is zero, CGET51 does nothing.
50 * It must be at least zero.
51 *
52 * A (input) COMPLEX array, dimension (LDA, N)
53 * The original (unfactored) matrix.
54 *
55 * LDA (input) INTEGER
56 * The leading dimension of A. It must be at least 1
57 * and at least N.
58 *
59 * B (input) COMPLEX array, dimension (LDB, N)
60 * The factored matrix.
61 *
62 * LDB (input) INTEGER
63 * The leading dimension of B. It must be at least 1
64 * and at least N.
65 *
66 * U (input) COMPLEX array, dimension (LDU, N)
67 * The unitary matrix on the left-hand side in the
68 * decomposition.
69 * Not referenced if ITYPE=2
70 *
71 * LDU (input) INTEGER
72 * The leading dimension of U. LDU must be at least N and
73 * at least 1.
74 *
75 * V (input) COMPLEX array, dimension (LDV, N)
76 * The unitary matrix on the left-hand side in the
77 * decomposition.
78 * Not referenced if ITYPE=2
79 *
80 * LDV (input) INTEGER
81 * The leading dimension of V. LDV must be at least N and
82 * at least 1.
83 *
84 * WORK (workspace) COMPLEX array, dimension (2*N**2)
85 *
86 * RWORK (workspace) REAL array, dimension (N)
87 *
88 * RESULT (output) REAL
89 * The values computed by the test specified by ITYPE. The
90 * value is currently limited to 1/ulp, to avoid overflow.
91 * Errors are flagged by RESULT=10/ulp.
92 *
93 * =====================================================================
94 *
95 * .. Parameters ..
96 REAL ZERO, ONE, TEN
97 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 10.0E+0 )
98 COMPLEX CZERO, CONE
99 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
100 $ CONE = ( 1.0E+0, 0.0E+0 ) )
101 * ..
102 * .. Local Scalars ..
103 INTEGER JCOL, JDIAG, JROW
104 REAL ANORM, ULP, UNFL, WNORM
105 * ..
106 * .. External Functions ..
107 REAL CLANGE, SLAMCH
108 EXTERNAL CLANGE, SLAMCH
109 * ..
110 * .. External Subroutines ..
111 EXTERNAL CGEMM, CLACPY
112 * ..
113 * .. Intrinsic Functions ..
114 INTRINSIC MAX, MIN, REAL
115 * ..
116 * .. Executable Statements ..
117 *
118 RESULT = ZERO
119 IF( N.LE.0 )
120 $ RETURN
121 *
122 * Constants
123 *
124 UNFL = SLAMCH( 'Safe minimum' )
125 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
126 *
127 * Some Error Checks
128 *
129 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
130 RESULT = TEN / ULP
131 RETURN
132 END IF
133 *
134 IF( ITYPE.LE.2 ) THEN
135 *
136 * Tests scaled by the norm(A)
137 *
138 ANORM = MAX( CLANGE( '1', N, N, A, LDA, RWORK ), UNFL )
139 *
140 IF( ITYPE.EQ.1 ) THEN
141 *
142 * ITYPE=1: Compute W = A - UBV'
143 *
144 CALL CLACPY( ' ', N, N, A, LDA, WORK, N )
145 CALL CGEMM( 'N', 'N', N, N, N, CONE, U, LDU, B, LDB, CZERO,
146 $ WORK( N**2+1 ), N )
147 *
148 CALL CGEMM( 'N', 'C', N, N, N, -CONE, WORK( N**2+1 ), N, V,
149 $ LDV, CONE, WORK, N )
150 *
151 ELSE
152 *
153 * ITYPE=2: Compute W = A - B
154 *
155 CALL CLACPY( ' ', N, N, B, LDB, WORK, N )
156 *
157 DO 20 JCOL = 1, N
158 DO 10 JROW = 1, N
159 WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) )
160 $ - A( JROW, JCOL )
161 10 CONTINUE
162 20 CONTINUE
163 END IF
164 *
165 * Compute norm(W)/ ( ulp*norm(A) )
166 *
167 WNORM = CLANGE( '1', N, N, WORK, N, RWORK )
168 *
169 IF( ANORM.GT.WNORM ) THEN
170 RESULT = ( WNORM / ANORM ) / ( N*ULP )
171 ELSE
172 IF( ANORM.LT.ONE ) THEN
173 RESULT = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
174 ELSE
175 RESULT = MIN( WNORM / ANORM, REAL( N ) ) / ( N*ULP )
176 END IF
177 END IF
178 *
179 ELSE
180 *
181 * Tests not scaled by norm(A)
182 *
183 * ITYPE=3: Compute UU' - I
184 *
185 CALL CGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO,
186 $ WORK, N )
187 *
188 DO 30 JDIAG = 1, N
189 WORK( ( N+1 )*( JDIAG-1 )+1 ) = WORK( ( N+1 )*( JDIAG-1 )+
190 $ 1 ) - CONE
191 30 CONTINUE
192 *
193 RESULT = MIN( CLANGE( '1', N, N, WORK, N, RWORK ),
194 $ REAL( N ) ) / ( N*ULP )
195 END IF
196 *
197 RETURN
198 *
199 * End of CGET51
200 *
201 END