1 SUBROUTINE CSPT03( UPLO, N, A, AINV, WORK, LDW, RWORK, RCOND,
2 $ RESID )
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 CHARACTER UPLO
10 INTEGER LDW, N
11 REAL RCOND, RESID
12 * ..
13 * .. Array Arguments ..
14 REAL RWORK( * )
15 COMPLEX A( * ), AINV( * ), WORK( LDW, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * CSPT03 computes the residual for a complex symmetric packed matrix
22 * times its inverse:
23 * norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
24 * where EPS is the machine epsilon.
25 *
26 * Arguments
27 * ==========
28 *
29 * UPLO (input) CHARACTER*1
30 * Specifies whether the upper or lower triangular part of the
31 * complex symmetric matrix A is stored:
32 * = 'U': Upper triangular
33 * = 'L': Lower triangular
34 *
35 * N (input) INTEGER
36 * The number of rows and columns of the matrix A. N >= 0.
37 *
38 * A (input) COMPLEX array, dimension (N*(N+1)/2)
39 * The original complex symmetric matrix A, stored as a packed
40 * triangular matrix.
41 *
42 * AINV (input) COMPLEX array, dimension (N*(N+1)/2)
43 * The (symmetric) inverse of the matrix A, stored as a packed
44 * triangular matrix.
45 *
46 * WORK (workspace) COMPLEX array, dimension (LDWORK,N)
47 *
48 * LDWORK (input) INTEGER
49 * The leading dimension of the array WORK. LDWORK >= max(1,N).
50 *
51 * RWORK (workspace) REAL array, dimension (N)
52 *
53 * RCOND (output) REAL
54 * The reciprocal of the condition number of A, computed as
55 * ( 1/norm(A) ) / norm(AINV).
56 *
57 * RESID (output) REAL
58 * norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 REAL ZERO, ONE
64 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
65 * ..
66 * .. Local Scalars ..
67 INTEGER I, ICOL, J, JCOL, K, KCOL, NALL
68 REAL AINVNM, ANORM, EPS
69 COMPLEX T
70 * ..
71 * .. External Functions ..
72 LOGICAL LSAME
73 REAL CLANGE, CLANSP, SLAMCH
74 COMPLEX CDOTU
75 EXTERNAL LSAME, CLANGE, CLANSP, SLAMCH, CDOTU
76 * ..
77 * .. Intrinsic Functions ..
78 INTRINSIC REAL
79 * ..
80 * .. Executable Statements ..
81 *
82 * Quick exit if N = 0.
83 *
84 IF( N.LE.0 ) THEN
85 RCOND = ONE
86 RESID = ZERO
87 RETURN
88 END IF
89 *
90 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
91 *
92 EPS = SLAMCH( 'Epsilon' )
93 ANORM = CLANSP( '1', UPLO, N, A, RWORK )
94 AINVNM = CLANSP( '1', UPLO, N, AINV, RWORK )
95 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
96 RCOND = ZERO
97 RESID = ONE / EPS
98 RETURN
99 END IF
100 RCOND = ( ONE/ANORM ) / AINVNM
101 *
102 * Case where both A and AINV are upper triangular:
103 * Each element of - A * AINV is computed by taking the dot product
104 * of a row of A with a column of AINV.
105 *
106 IF( LSAME( UPLO, 'U' ) ) THEN
107 DO 70 I = 1, N
108 ICOL = ( ( I-1 )*I ) / 2 + 1
109 *
110 * Code when J <= I
111 *
112 DO 30 J = 1, I
113 JCOL = ( ( J-1 )*J ) / 2 + 1
114 T = CDOTU( J, A( ICOL ), 1, AINV( JCOL ), 1 )
115 JCOL = JCOL + 2*J - 1
116 KCOL = ICOL - 1
117 DO 10 K = J + 1, I
118 T = T + A( KCOL+K )*AINV( JCOL )
119 JCOL = JCOL + K
120 10 CONTINUE
121 KCOL = KCOL + 2*I
122 DO 20 K = I + 1, N
123 T = T + A( KCOL )*AINV( JCOL )
124 KCOL = KCOL + K
125 JCOL = JCOL + K
126 20 CONTINUE
127 WORK( I, J ) = -T
128 30 CONTINUE
129 *
130 * Code when J > I
131 *
132 DO 60 J = I + 1, N
133 JCOL = ( ( J-1 )*J ) / 2 + 1
134 T = CDOTU( I, A( ICOL ), 1, AINV( JCOL ), 1 )
135 JCOL = JCOL - 1
136 KCOL = ICOL + 2*I - 1
137 DO 40 K = I + 1, J
138 T = T + A( KCOL )*AINV( JCOL+K )
139 KCOL = KCOL + K
140 40 CONTINUE
141 JCOL = JCOL + 2*J
142 DO 50 K = J + 1, N
143 T = T + A( KCOL )*AINV( JCOL )
144 KCOL = KCOL + K
145 JCOL = JCOL + K
146 50 CONTINUE
147 WORK( I, J ) = -T
148 60 CONTINUE
149 70 CONTINUE
150 ELSE
151 *
152 * Case where both A and AINV are lower triangular
153 *
154 NALL = ( N*( N+1 ) ) / 2
155 DO 140 I = 1, N
156 *
157 * Code when J <= I
158 *
159 ICOL = NALL - ( ( N-I+1 )*( N-I+2 ) ) / 2 + 1
160 DO 100 J = 1, I
161 JCOL = NALL - ( ( N-J )*( N-J+1 ) ) / 2 - ( N-I )
162 T = CDOTU( N-I+1, A( ICOL ), 1, AINV( JCOL ), 1 )
163 KCOL = I
164 JCOL = J
165 DO 80 K = 1, J - 1
166 T = T + A( KCOL )*AINV( JCOL )
167 JCOL = JCOL + N - K
168 KCOL = KCOL + N - K
169 80 CONTINUE
170 JCOL = JCOL - J
171 DO 90 K = J, I - 1
172 T = T + A( KCOL )*AINV( JCOL+K )
173 KCOL = KCOL + N - K
174 90 CONTINUE
175 WORK( I, J ) = -T
176 100 CONTINUE
177 *
178 * Code when J > I
179 *
180 ICOL = NALL - ( ( N-I )*( N-I+1 ) ) / 2
181 DO 130 J = I + 1, N
182 JCOL = NALL - ( ( N-J+1 )*( N-J+2 ) ) / 2 + 1
183 T = CDOTU( N-J+1, A( ICOL-N+J ), 1, AINV( JCOL ), 1 )
184 KCOL = I
185 JCOL = J
186 DO 110 K = 1, I - 1
187 T = T + A( KCOL )*AINV( JCOL )
188 JCOL = JCOL + N - K
189 KCOL = KCOL + N - K
190 110 CONTINUE
191 KCOL = KCOL - I
192 DO 120 K = I, J - 1
193 T = T + A( KCOL+K )*AINV( JCOL )
194 JCOL = JCOL + N - K
195 120 CONTINUE
196 WORK( I, J ) = -T
197 130 CONTINUE
198 140 CONTINUE
199 END IF
200 *
201 * Add the identity matrix to WORK .
202 *
203 DO 150 I = 1, N
204 WORK( I, I ) = WORK( I, I ) + ONE
205 150 CONTINUE
206 *
207 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
208 *
209 RESID = CLANGE( '1', N, N, WORK, LDW, RWORK )
210 *
211 RESID = ( ( RESID*RCOND )/EPS ) / REAL( N )
212 *
213 RETURN
214 *
215 * End of CSPT03
216 *
217 END
2 $ RESID )
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 CHARACTER UPLO
10 INTEGER LDW, N
11 REAL RCOND, RESID
12 * ..
13 * .. Array Arguments ..
14 REAL RWORK( * )
15 COMPLEX A( * ), AINV( * ), WORK( LDW, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * CSPT03 computes the residual for a complex symmetric packed matrix
22 * times its inverse:
23 * norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
24 * where EPS is the machine epsilon.
25 *
26 * Arguments
27 * ==========
28 *
29 * UPLO (input) CHARACTER*1
30 * Specifies whether the upper or lower triangular part of the
31 * complex symmetric matrix A is stored:
32 * = 'U': Upper triangular
33 * = 'L': Lower triangular
34 *
35 * N (input) INTEGER
36 * The number of rows and columns of the matrix A. N >= 0.
37 *
38 * A (input) COMPLEX array, dimension (N*(N+1)/2)
39 * The original complex symmetric matrix A, stored as a packed
40 * triangular matrix.
41 *
42 * AINV (input) COMPLEX array, dimension (N*(N+1)/2)
43 * The (symmetric) inverse of the matrix A, stored as a packed
44 * triangular matrix.
45 *
46 * WORK (workspace) COMPLEX array, dimension (LDWORK,N)
47 *
48 * LDWORK (input) INTEGER
49 * The leading dimension of the array WORK. LDWORK >= max(1,N).
50 *
51 * RWORK (workspace) REAL array, dimension (N)
52 *
53 * RCOND (output) REAL
54 * The reciprocal of the condition number of A, computed as
55 * ( 1/norm(A) ) / norm(AINV).
56 *
57 * RESID (output) REAL
58 * norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 REAL ZERO, ONE
64 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
65 * ..
66 * .. Local Scalars ..
67 INTEGER I, ICOL, J, JCOL, K, KCOL, NALL
68 REAL AINVNM, ANORM, EPS
69 COMPLEX T
70 * ..
71 * .. External Functions ..
72 LOGICAL LSAME
73 REAL CLANGE, CLANSP, SLAMCH
74 COMPLEX CDOTU
75 EXTERNAL LSAME, CLANGE, CLANSP, SLAMCH, CDOTU
76 * ..
77 * .. Intrinsic Functions ..
78 INTRINSIC REAL
79 * ..
80 * .. Executable Statements ..
81 *
82 * Quick exit if N = 0.
83 *
84 IF( N.LE.0 ) THEN
85 RCOND = ONE
86 RESID = ZERO
87 RETURN
88 END IF
89 *
90 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
91 *
92 EPS = SLAMCH( 'Epsilon' )
93 ANORM = CLANSP( '1', UPLO, N, A, RWORK )
94 AINVNM = CLANSP( '1', UPLO, N, AINV, RWORK )
95 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
96 RCOND = ZERO
97 RESID = ONE / EPS
98 RETURN
99 END IF
100 RCOND = ( ONE/ANORM ) / AINVNM
101 *
102 * Case where both A and AINV are upper triangular:
103 * Each element of - A * AINV is computed by taking the dot product
104 * of a row of A with a column of AINV.
105 *
106 IF( LSAME( UPLO, 'U' ) ) THEN
107 DO 70 I = 1, N
108 ICOL = ( ( I-1 )*I ) / 2 + 1
109 *
110 * Code when J <= I
111 *
112 DO 30 J = 1, I
113 JCOL = ( ( J-1 )*J ) / 2 + 1
114 T = CDOTU( J, A( ICOL ), 1, AINV( JCOL ), 1 )
115 JCOL = JCOL + 2*J - 1
116 KCOL = ICOL - 1
117 DO 10 K = J + 1, I
118 T = T + A( KCOL+K )*AINV( JCOL )
119 JCOL = JCOL + K
120 10 CONTINUE
121 KCOL = KCOL + 2*I
122 DO 20 K = I + 1, N
123 T = T + A( KCOL )*AINV( JCOL )
124 KCOL = KCOL + K
125 JCOL = JCOL + K
126 20 CONTINUE
127 WORK( I, J ) = -T
128 30 CONTINUE
129 *
130 * Code when J > I
131 *
132 DO 60 J = I + 1, N
133 JCOL = ( ( J-1 )*J ) / 2 + 1
134 T = CDOTU( I, A( ICOL ), 1, AINV( JCOL ), 1 )
135 JCOL = JCOL - 1
136 KCOL = ICOL + 2*I - 1
137 DO 40 K = I + 1, J
138 T = T + A( KCOL )*AINV( JCOL+K )
139 KCOL = KCOL + K
140 40 CONTINUE
141 JCOL = JCOL + 2*J
142 DO 50 K = J + 1, N
143 T = T + A( KCOL )*AINV( JCOL )
144 KCOL = KCOL + K
145 JCOL = JCOL + K
146 50 CONTINUE
147 WORK( I, J ) = -T
148 60 CONTINUE
149 70 CONTINUE
150 ELSE
151 *
152 * Case where both A and AINV are lower triangular
153 *
154 NALL = ( N*( N+1 ) ) / 2
155 DO 140 I = 1, N
156 *
157 * Code when J <= I
158 *
159 ICOL = NALL - ( ( N-I+1 )*( N-I+2 ) ) / 2 + 1
160 DO 100 J = 1, I
161 JCOL = NALL - ( ( N-J )*( N-J+1 ) ) / 2 - ( N-I )
162 T = CDOTU( N-I+1, A( ICOL ), 1, AINV( JCOL ), 1 )
163 KCOL = I
164 JCOL = J
165 DO 80 K = 1, J - 1
166 T = T + A( KCOL )*AINV( JCOL )
167 JCOL = JCOL + N - K
168 KCOL = KCOL + N - K
169 80 CONTINUE
170 JCOL = JCOL - J
171 DO 90 K = J, I - 1
172 T = T + A( KCOL )*AINV( JCOL+K )
173 KCOL = KCOL + N - K
174 90 CONTINUE
175 WORK( I, J ) = -T
176 100 CONTINUE
177 *
178 * Code when J > I
179 *
180 ICOL = NALL - ( ( N-I )*( N-I+1 ) ) / 2
181 DO 130 J = I + 1, N
182 JCOL = NALL - ( ( N-J+1 )*( N-J+2 ) ) / 2 + 1
183 T = CDOTU( N-J+1, A( ICOL-N+J ), 1, AINV( JCOL ), 1 )
184 KCOL = I
185 JCOL = J
186 DO 110 K = 1, I - 1
187 T = T + A( KCOL )*AINV( JCOL )
188 JCOL = JCOL + N - K
189 KCOL = KCOL + N - K
190 110 CONTINUE
191 KCOL = KCOL - I
192 DO 120 K = I, J - 1
193 T = T + A( KCOL+K )*AINV( JCOL )
194 JCOL = JCOL + N - K
195 120 CONTINUE
196 WORK( I, J ) = -T
197 130 CONTINUE
198 140 CONTINUE
199 END IF
200 *
201 * Add the identity matrix to WORK .
202 *
203 DO 150 I = 1, N
204 WORK( I, I ) = WORK( I, I ) + ONE
205 150 CONTINUE
206 *
207 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
208 *
209 RESID = CLANGE( '1', N, N, WORK, LDW, RWORK )
210 *
211 RESID = ( ( RESID*RCOND )/EPS ) / REAL( N )
212 *
213 RETURN
214 *
215 * End of CSPT03
216 *
217 END