 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ ctprfs()

 subroutine ctprfs ( character UPLO, character TRANS, character DIAG, integer N, integer NRHS, complex, dimension( * ) AP, complex, dimension( ldb, * ) B, integer LDB, complex, dimension( ldx, * ) X, integer LDX, real, dimension( * ) FERR, real, dimension( * ) BERR, complex, dimension( * ) WORK, real, dimension( * ) RWORK, integer INFO )

CTPRFS

Purpose:
``` CTPRFS provides error bounds and backward error estimates for the
solution to a system of linear equations with a triangular packed
coefficient matrix.

The solution matrix X must be computed by CTPTRS or some other
means before entering this routine.  CTPRFS does not do iterative
refinement because doing so cannot improve the backward error.```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 = 'U': A is upper triangular; = 'L': A is lower triangular.``` [in] TRANS ``` TRANS is CHARACTER*1 Specifies the form of the system of equations: = 'N': A * X = B (No transpose) = 'T': A**T * X = B (Transpose) = 'C': A**H * X = B (Conjugate transpose)``` [in] DIAG ``` DIAG is CHARACTER*1 = 'N': A is non-unit triangular; = 'U': A is unit triangular.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in] NRHS ``` NRHS is INTEGER The number of right hand sides, i.e., the number of columns of the matrices B and X. NRHS >= 0.``` [in] AP ``` AP is COMPLEX array, dimension (N*(N+1)/2) The upper or lower triangular matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows: if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. If DIAG = 'U', the diagonal elements of A are not referenced and are assumed to be 1.``` [in] B ``` B is COMPLEX array, dimension (LDB,NRHS) The right hand side matrix B.``` [in] LDB ``` LDB is INTEGER The leading dimension of the array B. LDB >= max(1,N).``` [in] X ``` X is COMPLEX array, dimension (LDX,NRHS) The solution matrix X.``` [in] LDX ``` LDX is INTEGER The leading dimension of the array X. LDX >= max(1,N).``` [out] FERR ``` FERR is REAL array, dimension (NRHS) The estimated forward error bound for each solution vector X(j) (the j-th column of the solution matrix X). If XTRUE is the true solution corresponding to X(j), FERR(j) is an estimated upper bound for the magnitude of the largest element in (X(j) - XTRUE) divided by the magnitude of the largest element in X(j). The estimate is as reliable as the estimate for RCOND, and is almost always a slight overestimate of the true error.``` [out] BERR ``` BERR is REAL array, dimension (NRHS) The componentwise relative backward error of each solution vector X(j) (i.e., the smallest relative change in any element of A or B that makes X(j) an exact solution).``` [out] WORK ` WORK is COMPLEX array, dimension (2*N)` [out] RWORK ` RWORK is REAL array, dimension (N)` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value```

Definition at line 172 of file ctprfs.f.

174*
175* -- LAPACK computational routine --
176* -- LAPACK is a software package provided by Univ. of Tennessee, --
177* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
178*
179* .. Scalar Arguments ..
180 CHARACTER DIAG, TRANS, UPLO
181 INTEGER INFO, LDB, LDX, N, NRHS
182* ..
183* .. Array Arguments ..
184 REAL BERR( * ), FERR( * ), RWORK( * )
185 COMPLEX AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
186* ..
187*
188* =====================================================================
189*
190* .. Parameters ..
191 REAL ZERO
192 parameter( zero = 0.0e+0 )
193 COMPLEX ONE
194 parameter( one = ( 1.0e+0, 0.0e+0 ) )
195* ..
196* .. Local Scalars ..
197 LOGICAL NOTRAN, NOUNIT, UPPER
198 CHARACTER TRANSN, TRANST
199 INTEGER I, J, K, KASE, KC, NZ
200 REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
201 COMPLEX ZDUM
202* ..
203* .. Local Arrays ..
204 INTEGER ISAVE( 3 )
205* ..
206* .. External Subroutines ..
207 EXTERNAL caxpy, ccopy, clacn2, ctpmv, ctpsv, 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 = 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 ) ) + nz*eps*rwork( i )
436 ELSE
437 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
438 \$ safe1
439 END IF
440 200 CONTINUE
441*
442 kase = 0
443 210 CONTINUE
444 CALL clacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
445 IF( kase.NE.0 ) THEN
446 IF( kase.EQ.1 ) THEN
447*
448* Multiply by diag(W)*inv(op(A)**H).
449*
450 CALL ctpsv( uplo, transt, diag, n, ap, work, 1 )
451 DO 220 i = 1, n
452 work( i ) = rwork( i )*work( i )
453 220 CONTINUE
454 ELSE
455*
456* Multiply by inv(op(A))*diag(W).
457*
458 DO 230 i = 1, n
459 work( i ) = rwork( i )*work( i )
460 230 CONTINUE
461 CALL ctpsv( uplo, transn, diag, n, ap, work, 1 )
462 END IF
463 GO TO 210
464 END IF
465*
466* Normalize error.
467*
468 lstres = zero
469 DO 240 i = 1, n
470 lstres = max( lstres, cabs1( x( i, j ) ) )
471 240 CONTINUE
472 IF( lstres.NE.zero )
473 \$ ferr( j ) = ferr( j ) / lstres
474*
475 250 CONTINUE
476*
477 RETURN
478*
479* End of CTPRFS
480*
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
subroutine ctpsv(UPLO, TRANS, DIAG, N, AP, X, INCX)
CTPSV
Definition: ctpsv.f:144
subroutine ctpmv(UPLO, TRANS, DIAG, N, AP, X, INCX)
CTPMV
Definition: ctpmv.f:142
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:133
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: