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