LAPACK 3.12.0
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*> \htmlonly
9*> Download ZTPRFS + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztprfs.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztprfs.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztprfs.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZTPRFS( UPLO, TRANS, DIAG, N, NRHS, AP, B, LDB, X, LDX,
22* FERR, BERR, WORK, RWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER DIAG, TRANS, UPLO
26* INTEGER INFO, LDB, LDX, N, NRHS
27* ..
28* .. Array Arguments ..
29* DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
30* COMPLEX*16 AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> ZTPRFS provides error bounds and backward error estimates for the
40*> solution to a system of linear equations with a triangular packed
41*> coefficient matrix.
42*>
43*> The solution matrix X must be computed by ZTPTRS or some other
44*> means before entering this routine. ZTPRFS does not do iterative
45*> refinement because doing so cannot improve the backward error.
46*> \endverbatim
47*
48* Arguments:
49* ==========
50*
51*> \param[in] UPLO
52*> \verbatim
53*> UPLO is CHARACTER*1
54*> = 'U': A is upper triangular;
55*> = 'L': A is lower triangular.
56*> \endverbatim
57*>
58*> \param[in] TRANS
59*> \verbatim
60*> TRANS is CHARACTER*1
61*> Specifies the form of the system of equations:
62*> = 'N': A * X = B (No transpose)
63*> = 'T': A**T * X = B (Transpose)
64*> = 'C': A**H * X = B (Conjugate transpose)
65*> \endverbatim
66*>
67*> \param[in] DIAG
68*> \verbatim
69*> DIAG is CHARACTER*1
70*> = 'N': A is non-unit triangular;
71*> = 'U': A is unit triangular.
72*> \endverbatim
73*>
74*> \param[in] N
75*> \verbatim
76*> N is INTEGER
77*> The order of the matrix A. N >= 0.
78*> \endverbatim
79*>
80*> \param[in] NRHS
81*> \verbatim
82*> NRHS is INTEGER
83*> The number of right hand sides, i.e., the number of columns
84*> of the matrices B and X. NRHS >= 0.
85*> \endverbatim
86*>
87*> \param[in] AP
88*> \verbatim
89*> AP is COMPLEX*16 array, dimension (N*(N+1)/2)
90*> The upper or lower triangular matrix A, packed columnwise in
91*> a linear array. The j-th column of A is stored in the array
92*> AP as follows:
93*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
94*> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
95*> If DIAG = 'U', the diagonal elements of A are not referenced
96*> and are assumed to be 1.
97*> \endverbatim
98*>
99*> \param[in] B
100*> \verbatim
101*> B is COMPLEX*16 array, dimension (LDB,NRHS)
102*> The right hand side matrix B.
103*> \endverbatim
104*>
105*> \param[in] LDB
106*> \verbatim
107*> LDB is INTEGER
108*> The leading dimension of the array B. LDB >= max(1,N).
109*> \endverbatim
110*>
111*> \param[in] X
112*> \verbatim
113*> X is COMPLEX*16 array, dimension (LDX,NRHS)
114*> The solution matrix X.
115*> \endverbatim
116*>
117*> \param[in] LDX
118*> \verbatim
119*> LDX is INTEGER
120*> The leading dimension of the array X. LDX >= max(1,N).
121*> \endverbatim
122*>
123*> \param[out] FERR
124*> \verbatim
125*> FERR is DOUBLE PRECISION array, dimension (NRHS)
126*> The estimated forward error bound for each solution vector
127*> X(j) (the j-th column of the solution matrix X).
128*> If XTRUE is the true solution corresponding to X(j), FERR(j)
129*> is an estimated upper bound for the magnitude of the largest
130*> element in (X(j) - XTRUE) divided by the magnitude of the
131*> largest element in X(j). The estimate is as reliable as
132*> the estimate for RCOND, and is almost always a slight
133*> overestimate of the true error.
134*> \endverbatim
135*>
136*> \param[out] BERR
137*> \verbatim
138*> BERR is DOUBLE PRECISION array, dimension (NRHS)
139*> The componentwise relative backward error of each solution
140*> vector X(j) (i.e., the smallest relative change in
141*> any element of A or B that makes X(j) an exact solution).
142*> \endverbatim
143*>
144*> \param[out] WORK
145*> \verbatim
146*> WORK is COMPLEX*16 array, dimension (2*N)
147*> \endverbatim
148*>
149*> \param[out] RWORK
150*> \verbatim
151*> RWORK is DOUBLE PRECISION array, dimension (N)
152*> \endverbatim
153*>
154*> \param[out] INFO
155*> \verbatim
156*> INFO is INTEGER
157*> = 0: successful exit
158*> < 0: if INFO = -i, the i-th argument had an illegal value
159*> \endverbatim
160*
161* Authors:
162* ========
163*
164*> \author Univ. of Tennessee
165*> \author Univ. of California Berkeley
166*> \author Univ. of Colorado Denver
167*> \author NAG Ltd.
168*
169*> \ingroup tprfs
170*
171* =====================================================================
172 SUBROUTINE ztprfs( UPLO, TRANS, DIAG, N, NRHS, AP, B, LDB, X, LDX,
173 $ FERR, BERR, WORK, RWORK, INFO )
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*
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:133
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:174
ztpsv
subroutine ztpsv(uplo, trans, diag, n, ap, x, incx)
ZTPSV
Definition ztpsv.f:144