LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ cspt03()

subroutine cspt03 ( character  uplo,
integer  n,
complex, dimension( * )  a,
complex, dimension( * )  ainv,
complex, dimension( ldw, * )  work,
integer  ldw,
real, dimension( * )  rwork,
real  rcond,
real  resid 
)

CSPT03

Purpose:
 CSPT03 computes the residual for a complex symmetric packed matrix
 times its inverse:
    norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
 where EPS is the machine epsilon.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          complex symmetric matrix A is stored:
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]N
          N is INTEGER
          The number of rows and columns of the matrix A.  N >= 0.
[in]A
          A is COMPLEX array, dimension (N*(N+1)/2)
          The original complex symmetric matrix A, stored as a packed
          triangular matrix.
[in]AINV
          AINV is COMPLEX array, dimension (N*(N+1)/2)
          The (symmetric) inverse of the matrix A, stored as a packed
          triangular matrix.
[out]WORK
          WORK is COMPLEX array, dimension (LDW,N)
[in]LDW
          LDW is INTEGER
          The leading dimension of the array WORK.  LDW >= max(1,N).
[out]RWORK
          RWORK is REAL array, dimension (N)
[out]RCOND
          RCOND is REAL
          The reciprocal of the condition number of A, computed as
          ( 1/norm(A) ) / norm(AINV).
[out]RESID
          RESID is REAL
          norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 108 of file cspt03.f.

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*
complex function cdotu(n, cx, incx, cy, incy)
CDOTU
Definition cdotu.f:83
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clange.f:115
real function clansp(norm, uplo, n, ap, work)
CLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clansp.f:115
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the caller graph for this function: