LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
ctrrfs.f
Go to the documentation of this file.
1*> \brief \b CTRRFS
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download CTRRFS + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctrrfs.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctrrfs.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctrrfs.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
20* LDX, FERR, BERR, WORK, RWORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER DIAG, TRANS, UPLO
24* INTEGER INFO, LDA, LDB, LDX, N, NRHS
25* ..
26* .. Array Arguments ..
27* REAL BERR( * ), FERR( * ), RWORK( * )
28* COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ),
29* $ X( LDX, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> CTRRFS provides error bounds and backward error estimates for the
39*> solution to a system of linear equations with a triangular
40*> coefficient matrix.
41*>
42*> The solution matrix X must be computed by CTRTRS or some other
43*> means before entering this routine. CTRRFS does not do iterative
44*> refinement because doing so cannot improve the backward error.
45*> \endverbatim
46*
47* Arguments:
48* ==========
49*
50*> \param[in] UPLO
51*> \verbatim
52*> UPLO is CHARACTER*1
53*> = 'U': A is upper triangular;
54*> = 'L': A is lower triangular.
55*> \endverbatim
56*>
57*> \param[in] TRANS
58*> \verbatim
59*> TRANS is CHARACTER*1
60*> Specifies the form of the system of equations:
61*> = 'N': A * X = B (No transpose)
62*> = 'T': A**T * X = B (Transpose)
63*> = 'C': A**H * X = B (Conjugate transpose)
64*> \endverbatim
65*>
66*> \param[in] DIAG
67*> \verbatim
68*> DIAG is CHARACTER*1
69*> = 'N': A is non-unit triangular;
70*> = 'U': A is unit triangular.
71*> \endverbatim
72*>
73*> \param[in] N
74*> \verbatim
75*> N is INTEGER
76*> The order of the matrix A. N >= 0.
77*> \endverbatim
78*>
79*> \param[in] NRHS
80*> \verbatim
81*> NRHS is INTEGER
82*> The number of right hand sides, i.e., the number of columns
83*> of the matrices B and X. NRHS >= 0.
84*> \endverbatim
85*>
86*> \param[in] A
87*> \verbatim
88*> A is COMPLEX array, dimension (LDA,N)
89*> The triangular matrix A. If UPLO = 'U', the leading N-by-N
90*> upper triangular part of the array A contains the upper
91*> triangular matrix, and the strictly lower triangular part of
92*> A is not referenced. If UPLO = 'L', the leading N-by-N lower
93*> triangular part of the array A contains the lower triangular
94*> matrix, and the strictly upper triangular part of A is not
95*> referenced. If DIAG = 'U', the diagonal elements of A are
96*> also not referenced and are assumed to be 1.
97*> \endverbatim
98*>
99*> \param[in] LDA
100*> \verbatim
101*> LDA is INTEGER
102*> The leading dimension of the array A. LDA >= max(1,N).
103*> \endverbatim
104*>
105*> \param[in] B
106*> \verbatim
107*> B is COMPLEX array, dimension (LDB,NRHS)
108*> The right hand side matrix B.
109*> \endverbatim
110*>
111*> \param[in] LDB
112*> \verbatim
113*> LDB is INTEGER
114*> The leading dimension of the array B. LDB >= max(1,N).
115*> \endverbatim
116*>
117*> \param[in] X
118*> \verbatim
119*> X is COMPLEX array, dimension (LDX,NRHS)
120*> The solution matrix X.
121*> \endverbatim
122*>
123*> \param[in] LDX
124*> \verbatim
125*> LDX is INTEGER
126*> The leading dimension of the array X. LDX >= max(1,N).
127*> \endverbatim
128*>
129*> \param[out] FERR
130*> \verbatim
131*> FERR is REAL array, dimension (NRHS)
132*> The estimated forward error bound for each solution vector
133*> X(j) (the j-th column of the solution matrix X).
134*> If XTRUE is the true solution corresponding to X(j), FERR(j)
135*> is an estimated upper bound for the magnitude of the largest
136*> element in (X(j) - XTRUE) divided by the magnitude of the
137*> largest element in X(j). The estimate is as reliable as
138*> the estimate for RCOND, and is almost always a slight
139*> overestimate of the true error.
140*> \endverbatim
141*>
142*> \param[out] BERR
143*> \verbatim
144*> BERR is REAL array, dimension (NRHS)
145*> The componentwise relative backward error of each solution
146*> vector X(j) (i.e., the smallest relative change in
147*> any element of A or B that makes X(j) an exact solution).
148*> \endverbatim
149*>
150*> \param[out] WORK
151*> \verbatim
152*> WORK is COMPLEX array, dimension (2*N)
153*> \endverbatim
154*>
155*> \param[out] RWORK
156*> \verbatim
157*> RWORK is REAL array, dimension (N)
158*> \endverbatim
159*>
160*> \param[out] INFO
161*> \verbatim
162*> INFO is INTEGER
163*> = 0: successful exit
164*> < 0: if INFO = -i, the i-th argument had an illegal value
165*> \endverbatim
166*
167* Authors:
168* ========
169*
170*> \author Univ. of Tennessee
171*> \author Univ. of California Berkeley
172*> \author Univ. of Colorado Denver
173*> \author NAG Ltd.
174*
175*> \ingroup trrfs
176*
177* =====================================================================
178 SUBROUTINE ctrrfs( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB,
179 $ X,
180 $ LDX, FERR, BERR, WORK, RWORK, INFO )
181*
182* -- LAPACK computational routine --
183* -- LAPACK is a software package provided by Univ. of Tennessee, --
184* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
185*
186* .. Scalar Arguments ..
187 CHARACTER DIAG, TRANS, UPLO
188 INTEGER INFO, LDA, LDB, LDX, N, NRHS
189* ..
190* .. Array Arguments ..
191 REAL BERR( * ), FERR( * ), RWORK( * )
192 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ),
193 $ x( ldx, * )
194* ..
195*
196* =====================================================================
197*
198* .. Parameters ..
199 REAL ZERO
200 PARAMETER ( ZERO = 0.0e+0 )
201 COMPLEX ONE
202 parameter( one = ( 1.0e+0, 0.0e+0 ) )
203* ..
204* .. Local Scalars ..
205 LOGICAL NOTRAN, NOUNIT, UPPER
206 CHARACTER TRANSN, TRANST
207 INTEGER I, J, K, KASE, NZ
208 REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
209 COMPLEX ZDUM
210* ..
211* .. Local Arrays ..
212 INTEGER ISAVE( 3 )
213* ..
214* .. External Subroutines ..
215 EXTERNAL caxpy, ccopy, clacn2, ctrmv, ctrsv,
216 $ xerbla
217* ..
218* .. Intrinsic Functions ..
219 INTRINSIC abs, aimag, max, real
220* ..
221* .. External Functions ..
222 LOGICAL LSAME
223 REAL SLAMCH
224 EXTERNAL lsame, slamch
225* ..
226* .. Statement Functions ..
227 REAL CABS1
228* ..
229* .. Statement Function definitions ..
230 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( 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( 'CTRRFS', -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 = slamch( 'Epsilon' )
286 safmin = slamch( 'Safe minimum' )
287 safe1 = real( 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 ccopy( n, x( 1, j ), 1, work, 1 )
298 CALL ctrmv( uplo, trans, diag, n, a, lda, work, 1 )
299 CALL caxpy( 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 CLACN2 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 ) ) + real( nz )*
431 $ eps*rwork( i )
432 ELSE
433 rwork( i ) = cabs1( work( i ) ) + real( nz )*
434 $ eps*rwork( i ) + safe1
435 END IF
436 200 CONTINUE
437*
438 kase = 0
439 210 CONTINUE
440 CALL clacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
441 IF( kase.NE.0 ) THEN
442 IF( kase.EQ.1 ) THEN
443*
444* Multiply by diag(W)*inv(op(A)**H).
445*
446 CALL ctrsv( uplo, transt, diag, n, a, lda, work, 1 )
447 DO 220 i = 1, n
448 work( i ) = rwork( i )*work( i )
449 220 CONTINUE
450 ELSE
451*
452* Multiply by inv(op(A))*diag(W).
453*
454 DO 230 i = 1, n
455 work( i ) = rwork( i )*work( i )
456 230 CONTINUE
457 CALL ctrsv( uplo, transn, diag, n, a, lda, work, 1 )
458 END IF
459 GO TO 210
460 END IF
461*
462* Normalize error.
463*
464 lstres = zero
465 DO 240 i = 1, n
466 lstres = max( lstres, cabs1( x( i, j ) ) )
467 240 CONTINUE
468 IF( lstres.NE.zero )
469 $ ferr( j ) = ferr( j ) / lstres
470*
471 250 CONTINUE
472*
473 RETURN
474*
475* End of CTRRFS
476*
477 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
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:131
subroutine ctrmv(uplo, trans, diag, n, a, lda, x, incx)
CTRMV
Definition ctrmv.f:147
subroutine ctrrfs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, x, ldx, ferr, berr, work, rwork, info)
CTRRFS
Definition ctrrfs.f:181
subroutine ctrsv(uplo, trans, diag, n, a, lda, x, incx)
CTRSV
Definition ctrsv.f:149