LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 complex16OTHERcomputational
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:88
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine ztrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
ZTRMV
Definition: ztrmv.f:147
subroutine ztrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
ZTRSV
Definition: ztrsv.f:149
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
subroutine ztrrfs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
ZTRRFS
Definition: ztrrfs.f:182