LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ztprfs.f
Go to the documentation of this file.
1*> \brief \b ZTPRFS
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download ZTPRFS + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztprfs.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztprfs.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztprfs.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE ZTPRFS( 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* DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
28* COMPLEX*16 AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> ZTPRFS 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 ZTPTRS or some other
42*> means before entering this routine. ZTPRFS 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*16 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*16 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*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 array, dimension (2*N)
145*> \endverbatim
146*>
147*> \param[out] RWORK
148*> \verbatim
149*> RWORK is DOUBLE PRECISION 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 ztprfs( 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 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
184 COMPLEX*16 AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
185* ..
186*
187* =====================================================================
188*
189* .. Parameters ..
190 DOUBLE PRECISION ZERO
191 PARAMETER ( ZERO = 0.0d+0 )
192 COMPLEX*16 ONE
193 parameter( one = ( 1.0d+0, 0.0d+0 ) )
194* ..
195* .. Local Scalars ..
196 LOGICAL NOTRAN, NOUNIT, UPPER
197 CHARACTER TRANSN, TRANST
198 INTEGER I, J, K, KASE, KC, NZ
199 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
200 COMPLEX*16 ZDUM
201* ..
202* .. Local Arrays ..
203 INTEGER ISAVE( 3 )
204* ..
205* .. External Subroutines ..
206 EXTERNAL xerbla, zaxpy, zcopy, zlacn2, ztpmv,
207 $ 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*
481 END
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:131
ztpmv
subroutine ztpmv(uplo, trans, diag, n, ap, x, incx)
ZTPMV
Definition ztpmv.f:142
ztprfs
subroutine ztprfs(uplo, trans, diag, n, nrhs, ap, b, ldb, x, ldx, ferr, berr, work, rwork, info)
ZTPRFS
Definition ztprfs.f:173
ztpsv
subroutine ztpsv(uplo, trans, diag, n, ap, x, incx)
ZTPSV
Definition ztpsv.f:144