LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cspt03.f
Go to the documentation of this file.
1*> \brief \b CSPT03
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE CSPT03( UPLO, N, A, AINV, WORK, LDW, RWORK, RCOND,
12* RESID )
13*
14* .. Scalar Arguments ..
15* CHARACTER UPLO
16* INTEGER LDW, N
17* REAL RCOND, RESID
18* ..
19* .. Array Arguments ..
20* REAL RWORK( * )
21* COMPLEX A( * ), AINV( * ), WORK( LDW, * )
22* ..
23*
24*
25*> \par Purpose:
26* =============
27*>
28*> \verbatim
29*>
30*> CSPT03 computes the residual for a complex symmetric packed matrix
31*> times its inverse:
32*> norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
33*> where EPS is the machine epsilon.
34*> \endverbatim
35*
36* Arguments:
37* ==========
38*
39*> \param[in] UPLO
40*> \verbatim
41*> UPLO is CHARACTER*1
42*> Specifies whether the upper or lower triangular part of the
43*> complex symmetric matrix A is stored:
44*> = 'U': Upper triangular
45*> = 'L': Lower triangular
46*> \endverbatim
47*>
48*> \param[in] N
49*> \verbatim
50*> N is INTEGER
51*> The number of rows and columns of the matrix A. N >= 0.
52*> \endverbatim
53*>
54*> \param[in] A
55*> \verbatim
56*> A is COMPLEX array, dimension (N*(N+1)/2)
57*> The original complex symmetric matrix A, stored as a packed
58*> triangular matrix.
59*> \endverbatim
60*>
61*> \param[in] AINV
62*> \verbatim
63*> AINV is COMPLEX array, dimension (N*(N+1)/2)
64*> The (symmetric) inverse of the matrix A, stored as a packed
65*> triangular matrix.
66*> \endverbatim
67*>
68*> \param[out] WORK
69*> \verbatim
70*> WORK is COMPLEX array, dimension (LDW,N)
71*> \endverbatim
72*>
73*> \param[in] LDW
74*> \verbatim
75*> LDW is INTEGER
76*> The leading dimension of the array WORK. LDW >= max(1,N).
77*> \endverbatim
78*>
79*> \param[out] RWORK
80*> \verbatim
81*> RWORK is REAL array, dimension (N)
82*> \endverbatim
83*>
84*> \param[out] RCOND
85*> \verbatim
86*> RCOND is REAL
87*> The reciprocal of the condition number of A, computed as
88*> ( 1/norm(A) ) / norm(AINV).
89*> \endverbatim
90*>
91*> \param[out] RESID
92*> \verbatim
93*> RESID is REAL
94*> norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
95*> \endverbatim
96*
97* Authors:
98* ========
99*
100*> \author Univ. of Tennessee
101*> \author Univ. of California Berkeley
102*> \author Univ. of Colorado Denver
103*> \author NAG Ltd.
104*
105*> \ingroup complex_lin
106*
107* =====================================================================
108 SUBROUTINE cspt03( UPLO, N, A, AINV, WORK, LDW, RWORK, RCOND,
109 $ RESID )
110*
111* -- LAPACK test routine --
112* -- LAPACK is a software package provided by Univ. of Tennessee, --
113* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
114*
115* .. Scalar Arguments ..
116 CHARACTER UPLO
117 INTEGER LDW, N
118 REAL RCOND, RESID
119* ..
120* .. Array Arguments ..
121 REAL RWORK( * )
122 COMPLEX A( * ), AINV( * ), WORK( LDW, * )
123* ..
124*
125* =====================================================================
126*
127* .. Parameters ..
128 REAL ZERO, ONE
129 parameter( zero = 0.0e+0, one = 1.0e+0 )
130* ..
131* .. Local Scalars ..
132 INTEGER I, ICOL, J, JCOL, K, KCOL, NALL
133 REAL AINVNM, ANORM, EPS
134 COMPLEX T
135* ..
136* .. External Functions ..
137 LOGICAL LSAME
138 REAL CLANGE, CLANSP, SLAMCH
139 COMPLEX CDOTU
140 EXTERNAL lsame, clange, clansp, slamch, cdotu
141* ..
142* .. Intrinsic Functions ..
143 INTRINSIC real
144* ..
145* .. Executable Statements ..
146*
147* Quick exit if N = 0.
148*
149 IF( n.LE.0 ) THEN
150 rcond = one
151 resid = zero
152 RETURN
153 END IF
154*
155* Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
156*
157 eps = slamch( 'Epsilon' )
158 anorm = clansp( '1', uplo, n, a, rwork )
159 ainvnm = clansp( '1', uplo, n, ainv, rwork )
160 IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
161 rcond = zero
162 resid = one / eps
163 RETURN
164 END IF
165 rcond = ( one/anorm ) / ainvnm
166*
167* Case where both A and AINV are upper triangular:
168* Each element of - A * AINV is computed by taking the dot product
169* of a row of A with a column of AINV.
170*
171 IF( lsame( uplo, 'U' ) ) THEN
172 DO 70 i = 1, n
173 icol = ( ( i-1 )*i ) / 2 + 1
174*
175* Code when J <= I
176*
177 DO 30 j = 1, i
178 jcol = ( ( j-1 )*j ) / 2 + 1
179 t = cdotu( j, a( icol ), 1, ainv( jcol ), 1 )
180 jcol = jcol + 2*j - 1
181 kcol = icol - 1
182 DO 10 k = j + 1, i
183 t = t + a( kcol+k )*ainv( jcol )
184 jcol = jcol + k
185 10 CONTINUE
186 kcol = kcol + 2*i
187 DO 20 k = i + 1, n
188 t = t + a( kcol )*ainv( jcol )
189 kcol = kcol + k
190 jcol = jcol + k
191 20 CONTINUE
192 work( i, j ) = -t
193 30 CONTINUE
194*
195* Code when J > I
196*
197 DO 60 j = i + 1, n
198 jcol = ( ( j-1 )*j ) / 2 + 1
199 t = cdotu( i, a( icol ), 1, ainv( jcol ), 1 )
200 jcol = jcol - 1
201 kcol = icol + 2*i - 1
202 DO 40 k = i + 1, j
203 t = t + a( kcol )*ainv( jcol+k )
204 kcol = kcol + k
205 40 CONTINUE
206 jcol = jcol + 2*j
207 DO 50 k = j + 1, n
208 t = t + a( kcol )*ainv( jcol )
209 kcol = kcol + k
210 jcol = jcol + k
211 50 CONTINUE
212 work( i, j ) = -t
213 60 CONTINUE
214 70 CONTINUE
215 ELSE
216*
217* Case where both A and AINV are lower triangular
218*
219 nall = ( n*( n+1 ) ) / 2
220 DO 140 i = 1, n
221*
222* Code when J <= I
223*
224 icol = nall - ( ( n-i+1 )*( n-i+2 ) ) / 2 + 1
225 DO 100 j = 1, i
226 jcol = nall - ( ( n-j )*( n-j+1 ) ) / 2 - ( n-i )
227 t = cdotu( n-i+1, a( icol ), 1, ainv( jcol ), 1 )
228 kcol = i
229 jcol = j
230 DO 80 k = 1, j - 1
231 t = t + a( kcol )*ainv( jcol )
232 jcol = jcol + n - k
233 kcol = kcol + n - k
234 80 CONTINUE
235 jcol = jcol - j
236 DO 90 k = j, i - 1
237 t = t + a( kcol )*ainv( jcol+k )
238 kcol = kcol + n - k
239 90 CONTINUE
240 work( i, j ) = -t
241 100 CONTINUE
242*
243* Code when J > I
244*
245 icol = nall - ( ( n-i )*( n-i+1 ) ) / 2
246 DO 130 j = i + 1, n
247 jcol = nall - ( ( n-j+1 )*( n-j+2 ) ) / 2 + 1
248 t = cdotu( n-j+1, a( icol-n+j ), 1, ainv( jcol ), 1 )
249 kcol = i
250 jcol = j
251 DO 110 k = 1, i - 1
252 t = t + a( kcol )*ainv( jcol )
253 jcol = jcol + n - k
254 kcol = kcol + n - k
255 110 CONTINUE
256 kcol = kcol - i
257 DO 120 k = i, j - 1
258 t = t + a( kcol+k )*ainv( jcol )
259 jcol = jcol + n - k
260 120 CONTINUE
261 work( i, j ) = -t
262 130 CONTINUE
263 140 CONTINUE
264 END IF
265*
266* Add the identity matrix to WORK .
267*
268 DO 150 i = 1, n
269 work( i, i ) = work( i, i ) + one
270 150 CONTINUE
271*
272* Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
273*
274 resid = clange( '1', n, n, work, ldw, rwork )
275*
276 resid = ( ( resid*rcond )/eps ) / real( n )
277*
278 RETURN
279*
280* End of CSPT03
281*
282 END
subroutine cspt03(uplo, n, a, ainv, work, ldw, rwork, rcond, resid)
CSPT03
Definition cspt03.f:110