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

◆ ztprfs()

subroutine ztprfs ( character  uplo,
character  trans,
character  diag,
integer  n,
integer  nrhs,
complex*16, dimension( * )  ap,
complex*16, dimension( ldb, * )  b,
integer  ldb,
complex*16, dimension( ldx, * )  x,
integer  ldx,
double precision, dimension( * )  ferr,
double precision, dimension( * )  berr,
complex*16, dimension( * )  work,
double precision, dimension( * )  rwork,
integer  info 
)

ZTPRFS

Download ZTPRFS + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZTPRFS 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 ZTPTRS or some other
 means before entering this routine.  ZTPRFS 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*16 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*16 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*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 172 of file ztprfs.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 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
185 COMPLEX*16 AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
186* ..
187*
188* =====================================================================
189*
190* .. Parameters ..
191 DOUBLE PRECISION ZERO
192 parameter( zero = 0.0d+0 )
193 COMPLEX*16 ONE
194 parameter( one = ( 1.0d+0, 0.0d+0 ) )
195* ..
196* .. Local Scalars ..
197 LOGICAL NOTRAN, NOUNIT, UPPER
198 CHARACTER TRANSN, TRANST
199 INTEGER I, J, K, KASE, KC, NZ
200 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
201 COMPLEX*16 ZDUM
202* ..
203* .. Local Arrays ..
204 INTEGER ISAVE( 3 )
205* ..
206* .. External Subroutines ..
207 EXTERNAL xerbla, zaxpy, zcopy, zlacn2, ztpmv, ztpsv
208* ..
209* .. Intrinsic Functions ..
210 INTRINSIC abs, dble, dimag, max
211* ..
212* .. External Functions ..
213 LOGICAL LSAME
214 DOUBLE PRECISION DLAMCH
215 EXTERNAL lsame, dlamch
216* ..
217* .. Statement Functions ..
218 DOUBLE PRECISION CABS1
219* ..
220* .. Statement Function definitions ..
221 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( 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( 'ZTPRFS', -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 = dlamch( 'Epsilon' )
275 safmin = dlamch( '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 zcopy( n, x( 1, j ), 1, work, 1 )
287 CALL ztpmv( uplo, trans, diag, n, ap, work, 1 )
288 CALL zaxpy( 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 ZLACN2 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 zlacn2( 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 ztpsv( 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 ztpsv( 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 ZTPRFS
480*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zaxpy
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zlacn2
subroutine zlacn2(n, v, x, est, kase, isave)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition zlacn2.f:133
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
ztpmv
subroutine ztpmv(uplo, trans, diag, n, ap, x, incx)
ZTPMV
Definition ztpmv.f:142
ztpsv
subroutine ztpsv(uplo, trans, diag, n, ap, x, incx)
ZTPSV
Definition ztpsv.f:144
Here is the call graph for this function:
Here is the caller graph for this function: