LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ctprfs.f
Go to the documentation of this file.
1*> \brief \b CTPRFS
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download CTPRFS + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctprfs.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctprfs.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctprfs.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CTPRFS( UPLO, TRANS, DIAG, N, NRHS, AP, B, LDB, X, LDX,
20* FERR, BERR, WORK, RWORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER DIAG, TRANS, UPLO
24* INTEGER INFO, LDB, LDX, N, NRHS
25* ..
26* .. Array Arguments ..
27* REAL BERR( * ), FERR( * ), RWORK( * )
28* COMPLEX AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> CTPRFS provides error bounds and backward error estimates for the
38*> solution to a system of linear equations with a triangular packed
39*> coefficient matrix.
40*>
41*> The solution matrix X must be computed by CTPTRS or some other
42*> means before entering this routine. CTPRFS does not do iterative
43*> refinement because doing so cannot improve the backward error.
44*> \endverbatim
45*
46* Arguments:
47* ==========
48*
49*> \param[in] UPLO
50*> \verbatim
51*> UPLO is CHARACTER*1
52*> = 'U': A is upper triangular;
53*> = 'L': A is lower triangular.
54*> \endverbatim
55*>
56*> \param[in] TRANS
57*> \verbatim
58*> TRANS is CHARACTER*1
59*> Specifies the form of the system of equations:
60*> = 'N': A * X = B (No transpose)
61*> = 'T': A**T * X = B (Transpose)
62*> = 'C': A**H * X = B (Conjugate transpose)
63*> \endverbatim
64*>
65*> \param[in] DIAG
66*> \verbatim
67*> DIAG is CHARACTER*1
68*> = 'N': A is non-unit triangular;
69*> = 'U': A is unit triangular.
70*> \endverbatim
71*>
72*> \param[in] N
73*> \verbatim
74*> N is INTEGER
75*> The order of the matrix A. N >= 0.
76*> \endverbatim
77*>
78*> \param[in] NRHS
79*> \verbatim
80*> NRHS is INTEGER
81*> The number of right hand sides, i.e., the number of columns
82*> of the matrices B and X. NRHS >= 0.
83*> \endverbatim
84*>
85*> \param[in] AP
86*> \verbatim
87*> AP is COMPLEX array, dimension (N*(N+1)/2)
88*> The upper or lower triangular matrix A, packed columnwise in
89*> a linear array. The j-th column of A is stored in the array
90*> AP as follows:
91*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
92*> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
93*> If DIAG = 'U', the diagonal elements of A are not referenced
94*> and are assumed to be 1.
95*> \endverbatim
96*>
97*> \param[in] B
98*> \verbatim
99*> B is COMPLEX array, dimension (LDB,NRHS)
100*> The right hand side matrix B.
101*> \endverbatim
102*>
103*> \param[in] LDB
104*> \verbatim
105*> LDB is INTEGER
106*> The leading dimension of the array B. LDB >= max(1,N).
107*> \endverbatim
108*>
109*> \param[in] X
110*> \verbatim
111*> X is COMPLEX array, dimension (LDX,NRHS)
112*> The solution matrix X.
113*> \endverbatim
114*>
115*> \param[in] LDX
116*> \verbatim
117*> LDX is INTEGER
118*> The leading dimension of the array X. LDX >= max(1,N).
119*> \endverbatim
120*>
121*> \param[out] FERR
122*> \verbatim
123*> FERR is REAL array, dimension (NRHS)
124*> The estimated forward error bound for each solution vector
125*> X(j) (the j-th column of the solution matrix X).
126*> If XTRUE is the true solution corresponding to X(j), FERR(j)
127*> is an estimated upper bound for the magnitude of the largest
128*> element in (X(j) - XTRUE) divided by the magnitude of the
129*> largest element in X(j). The estimate is as reliable as
130*> the estimate for RCOND, and is almost always a slight
131*> overestimate of the true error.
132*> \endverbatim
133*>
134*> \param[out] BERR
135*> \verbatim
136*> BERR is REAL array, dimension (NRHS)
137*> The componentwise relative backward error of each solution
138*> vector X(j) (i.e., the smallest relative change in
139*> any element of A or B that makes X(j) an exact solution).
140*> \endverbatim
141*>
142*> \param[out] WORK
143*> \verbatim
144*> WORK is COMPLEX array, dimension (2*N)
145*> \endverbatim
146*>
147*> \param[out] RWORK
148*> \verbatim
149*> RWORK is REAL array, dimension (N)
150*> \endverbatim
151*>
152*> \param[out] INFO
153*> \verbatim
154*> INFO is INTEGER
155*> = 0: successful exit
156*> < 0: if INFO = -i, the i-th argument had an illegal value
157*> \endverbatim
158*
159* Authors:
160* ========
161*
162*> \author Univ. of Tennessee
163*> \author Univ. of California Berkeley
164*> \author Univ. of Colorado Denver
165*> \author NAG Ltd.
166*
167*> \ingroup tprfs
168*
169* =====================================================================
170 SUBROUTINE ctprfs( UPLO, TRANS, DIAG, N, NRHS, AP, B, LDB, X,
171 $ LDX,
172 $ FERR, BERR, WORK, RWORK, INFO )
173*
174* -- LAPACK computational routine --
175* -- LAPACK is a software package provided by Univ. of Tennessee, --
176* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
177*
178* .. Scalar Arguments ..
179 CHARACTER DIAG, TRANS, UPLO
180 INTEGER INFO, LDB, LDX, N, NRHS
181* ..
182* .. Array Arguments ..
183 REAL BERR( * ), FERR( * ), RWORK( * )
184 COMPLEX AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
185* ..
186*
187* =====================================================================
188*
189* .. Parameters ..
190 REAL ZERO
191 PARAMETER ( ZERO = 0.0e+0 )
192 COMPLEX ONE
193 parameter( one = ( 1.0e+0, 0.0e+0 ) )
194* ..
195* .. Local Scalars ..
196 LOGICAL NOTRAN, NOUNIT, UPPER
197 CHARACTER TRANSN, TRANST
198 INTEGER I, J, K, KASE, KC, NZ
199 REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
200 COMPLEX ZDUM
201* ..
202* .. Local Arrays ..
203 INTEGER ISAVE( 3 )
204* ..
205* .. External Subroutines ..
206 EXTERNAL caxpy, ccopy, clacn2, ctpmv, ctpsv,
207 $ xerbla
208* ..
209* .. Intrinsic Functions ..
210 INTRINSIC abs, aimag, max, real
211* ..
212* .. External Functions ..
213 LOGICAL LSAME
214 REAL SLAMCH
215 EXTERNAL lsame, slamch
216* ..
217* .. Statement Functions ..
218 REAL CABS1
219* ..
220* .. Statement Function definitions ..
221 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
222* ..
223* .. Executable Statements ..
224*
225* Test the input parameters.
226*
227 info = 0
228 upper = lsame( uplo, 'U' )
229 notran = lsame( trans, 'N' )
230 nounit = lsame( diag, 'N' )
231*
232 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
233 info = -1
234 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
235 $ lsame( trans, 'C' ) ) THEN
236 info = -2
237 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
238 info = -3
239 ELSE IF( n.LT.0 ) THEN
240 info = -4
241 ELSE IF( nrhs.LT.0 ) THEN
242 info = -5
243 ELSE IF( ldb.LT.max( 1, n ) ) THEN
244 info = -8
245 ELSE IF( ldx.LT.max( 1, n ) ) THEN
246 info = -10
247 END IF
248 IF( info.NE.0 ) THEN
249 CALL xerbla( 'CTPRFS', -info )
250 RETURN
251 END IF
252*
253* Quick return if possible
254*
255 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
256 DO 10 j = 1, nrhs
257 ferr( j ) = zero
258 berr( j ) = zero
259 10 CONTINUE
260 RETURN
261 END IF
262*
263 IF( notran ) THEN
264 transn = 'N'
265 transt = 'C'
266 ELSE
267 transn = 'C'
268 transt = 'N'
269 END IF
270*
271* NZ = maximum number of nonzero elements in each row of A, plus 1
272*
273 nz = n + 1
274 eps = slamch( 'Epsilon' )
275 safmin = slamch( 'Safe minimum' )
276 safe1 = real( nz )*safmin
277 safe2 = safe1 / eps
278*
279* Do for each right hand side
280*
281 DO 250 j = 1, nrhs
282*
283* Compute residual R = B - op(A) * X,
284* where op(A) = A, A**T, or A**H, depending on TRANS.
285*
286 CALL ccopy( n, x( 1, j ), 1, work, 1 )
287 CALL ctpmv( uplo, trans, diag, n, ap, work, 1 )
288 CALL caxpy( n, -one, b( 1, j ), 1, work, 1 )
289*
290* Compute componentwise relative backward error from formula
291*
292* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
293*
294* where abs(Z) is the componentwise absolute value of the matrix
295* or vector Z. If the i-th component of the denominator is less
296* than SAFE2, then SAFE1 is added to the i-th components of the
297* numerator and denominator before dividing.
298*
299 DO 20 i = 1, n
300 rwork( i ) = cabs1( b( i, j ) )
301 20 CONTINUE
302*
303 IF( notran ) THEN
304*
305* Compute abs(A)*abs(X) + abs(B).
306*
307 IF( upper ) THEN
308 kc = 1
309 IF( nounit ) THEN
310 DO 40 k = 1, n
311 xk = cabs1( x( k, j ) )
312 DO 30 i = 1, k
313 rwork( i ) = rwork( i ) +
314 $ cabs1( ap( kc+i-1 ) )*xk
315 30 CONTINUE
316 kc = kc + k
317 40 CONTINUE
318 ELSE
319 DO 60 k = 1, n
320 xk = cabs1( x( k, j ) )
321 DO 50 i = 1, k - 1
322 rwork( i ) = rwork( i ) +
323 $ cabs1( ap( kc+i-1 ) )*xk
324 50 CONTINUE
325 rwork( k ) = rwork( k ) + xk
326 kc = kc + k
327 60 CONTINUE
328 END IF
329 ELSE
330 kc = 1
331 IF( nounit ) THEN
332 DO 80 k = 1, n
333 xk = cabs1( x( k, j ) )
334 DO 70 i = k, n
335 rwork( i ) = rwork( i ) +
336 $ cabs1( ap( kc+i-k ) )*xk
337 70 CONTINUE
338 kc = kc + n - k + 1
339 80 CONTINUE
340 ELSE
341 DO 100 k = 1, n
342 xk = cabs1( x( k, j ) )
343 DO 90 i = k + 1, n
344 rwork( i ) = rwork( i ) +
345 $ cabs1( ap( kc+i-k ) )*xk
346 90 CONTINUE
347 rwork( k ) = rwork( k ) + xk
348 kc = kc + n - k + 1
349 100 CONTINUE
350 END IF
351 END IF
352 ELSE
353*
354* Compute abs(A**H)*abs(X) + abs(B).
355*
356 IF( upper ) THEN
357 kc = 1
358 IF( nounit ) THEN
359 DO 120 k = 1, n
360 s = zero
361 DO 110 i = 1, k
362 s = s + cabs1( ap( kc+i-1 ) )*cabs1( x( i, j ) )
363 110 CONTINUE
364 rwork( k ) = rwork( k ) + s
365 kc = kc + k
366 120 CONTINUE
367 ELSE
368 DO 140 k = 1, n
369 s = cabs1( x( k, j ) )
370 DO 130 i = 1, k - 1
371 s = s + cabs1( ap( kc+i-1 ) )*cabs1( x( i, j ) )
372 130 CONTINUE
373 rwork( k ) = rwork( k ) + s
374 kc = kc + k
375 140 CONTINUE
376 END IF
377 ELSE
378 kc = 1
379 IF( nounit ) THEN
380 DO 160 k = 1, n
381 s = zero
382 DO 150 i = k, n
383 s = s + cabs1( ap( kc+i-k ) )*cabs1( x( i, j ) )
384 150 CONTINUE
385 rwork( k ) = rwork( k ) + s
386 kc = kc + n - k + 1
387 160 CONTINUE
388 ELSE
389 DO 180 k = 1, n
390 s = cabs1( x( k, j ) )
391 DO 170 i = k + 1, n
392 s = s + cabs1( ap( kc+i-k ) )*cabs1( x( i, j ) )
393 170 CONTINUE
394 rwork( k ) = rwork( k ) + s
395 kc = kc + n - k + 1
396 180 CONTINUE
397 END IF
398 END IF
399 END IF
400 s = zero
401 DO 190 i = 1, n
402 IF( rwork( i ).GT.safe2 ) THEN
403 s = max( s, cabs1( work( i ) ) / rwork( i ) )
404 ELSE
405 s = max( s, ( cabs1( work( i ) )+safe1 ) /
406 $ ( rwork( i )+safe1 ) )
407 END IF
408 190 CONTINUE
409 berr( j ) = s
410*
411* Bound error from formula
412*
413* norm(X - XTRUE) / norm(X) .le. FERR =
414* norm( abs(inv(op(A)))*
415* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
416*
417* where
418* norm(Z) is the magnitude of the largest component of Z
419* inv(op(A)) is the inverse of op(A)
420* abs(Z) is the componentwise absolute value of the matrix or
421* vector Z
422* NZ is the maximum number of nonzeros in any row of A, plus 1
423* EPS is machine epsilon
424*
425* The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
426* is incremented by SAFE1 if the i-th component of
427* abs(op(A))*abs(X) + abs(B) is less than SAFE2.
428*
429* Use CLACN2 to estimate the infinity-norm of the matrix
430* inv(op(A)) * diag(W),
431* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
432*
433 DO 200 i = 1, n
434 IF( rwork( i ).GT.safe2 ) THEN
435 rwork( i ) = cabs1( work( i ) ) + real( nz )*
436 $ eps*rwork( i )
437 ELSE
438 rwork( i ) = cabs1( work( i ) ) + real( nz )*
439 $ eps*rwork( i ) + safe1
440 END IF
441 200 CONTINUE
442*
443 kase = 0
444 210 CONTINUE
445 CALL clacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
446 IF( kase.NE.0 ) THEN
447 IF( kase.EQ.1 ) THEN
448*
449* Multiply by diag(W)*inv(op(A)**H).
450*
451 CALL ctpsv( uplo, transt, diag, n, ap, work, 1 )
452 DO 220 i = 1, n
453 work( i ) = rwork( i )*work( i )
454 220 CONTINUE
455 ELSE
456*
457* Multiply by inv(op(A))*diag(W).
458*
459 DO 230 i = 1, n
460 work( i ) = rwork( i )*work( i )
461 230 CONTINUE
462 CALL ctpsv( uplo, transn, diag, n, ap, work, 1 )
463 END IF
464 GO TO 210
465 END IF
466*
467* Normalize error.
468*
469 lstres = zero
470 DO 240 i = 1, n
471 lstres = max( lstres, cabs1( x( i, j ) ) )
472 240 CONTINUE
473 IF( lstres.NE.zero )
474 $ ferr( j ) = ferr( j ) / lstres
475*
476 250 CONTINUE
477*
478 RETURN
479*
480* End of CTPRFS
481*
482 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine clacn2(n, v, x, est, kase, isave)
CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition clacn2.f:131
subroutine ctpmv(uplo, trans, diag, n, ap, x, incx)
CTPMV
Definition ctpmv.f:142
subroutine ctprfs(uplo, trans, diag, n, nrhs, ap, b, ldb, x, ldx, ferr, berr, work, rwork, info)
CTPRFS
Definition ctprfs.f:173
subroutine ctpsv(uplo, trans, diag, n, ap, x, incx)
CTPSV
Definition ctpsv.f:144