LAPACK  3.6.1
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 *> \date November 2011
178 *
179 *> \ingroup complex16OTHERcomputational
180 *
181 * =====================================================================
182  SUBROUTINE ztrrfs( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
183  $ ldx, ferr, berr, work, rwork, info )
184 *
185 * -- LAPACK computational routine (version 3.4.0) --
186 * -- LAPACK is a software package provided by Univ. of Tennessee, --
187 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188 * November 2011
189 *
190 * .. Scalar Arguments ..
191  CHARACTER DIAG, TRANS, UPLO
192  INTEGER INFO, LDA, LDB, LDX, N, NRHS
193 * ..
194 * .. Array Arguments ..
195  DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
196  COMPLEX*16 A( lda, * ), B( ldb, * ), WORK( * ),
197  $ x( ldx, * )
198 * ..
199 *
200 * =====================================================================
201 *
202 * .. Parameters ..
203  DOUBLE PRECISION ZERO
204  parameter ( zero = 0.0d+0 )
205  COMPLEX*16 ONE
206  parameter ( one = ( 1.0d+0, 0.0d+0 ) )
207 * ..
208 * .. Local Scalars ..
209  LOGICAL NOTRAN, NOUNIT, UPPER
210  CHARACTER TRANSN, TRANST
211  INTEGER I, J, K, KASE, NZ
212  DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
213  COMPLEX*16 ZDUM
214 * ..
215 * .. Local Arrays ..
216  INTEGER ISAVE( 3 )
217 * ..
218 * .. External Subroutines ..
219  EXTERNAL xerbla, zaxpy, zcopy, zlacn2, ztrmv, ztrsv
220 * ..
221 * .. Intrinsic Functions ..
222  INTRINSIC abs, dble, dimag, max
223 * ..
224 * .. External Functions ..
225  LOGICAL LSAME
226  DOUBLE PRECISION DLAMCH
227  EXTERNAL lsame, dlamch
228 * ..
229 * .. Statement Functions ..
230  DOUBLE PRECISION CABS1
231 * ..
232 * .. Statement Function definitions ..
233  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
234 * ..
235 * .. Executable Statements ..
236 *
237 * Test the input parameters.
238 *
239  info = 0
240  upper = lsame( uplo, 'U' )
241  notran = lsame( trans, 'N' )
242  nounit = lsame( diag, 'N' )
243 *
244  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
245  info = -1
246  ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
247  $ lsame( trans, 'C' ) ) THEN
248  info = -2
249  ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
250  info = -3
251  ELSE IF( n.LT.0 ) THEN
252  info = -4
253  ELSE IF( nrhs.LT.0 ) THEN
254  info = -5
255  ELSE IF( lda.LT.max( 1, n ) ) THEN
256  info = -7
257  ELSE IF( ldb.LT.max( 1, n ) ) THEN
258  info = -9
259  ELSE IF( ldx.LT.max( 1, n ) ) THEN
260  info = -11
261  END IF
262  IF( info.NE.0 ) THEN
263  CALL xerbla( 'ZTRRFS', -info )
264  RETURN
265  END IF
266 *
267 * Quick return if possible
268 *
269  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
270  DO 10 j = 1, nrhs
271  ferr( j ) = zero
272  berr( j ) = zero
273  10 CONTINUE
274  RETURN
275  END IF
276 *
277  IF( notran ) THEN
278  transn = 'N'
279  transt = 'C'
280  ELSE
281  transn = 'C'
282  transt = 'N'
283  END IF
284 *
285 * NZ = maximum number of nonzero elements in each row of A, plus 1
286 *
287  nz = n + 1
288  eps = dlamch( 'Epsilon' )
289  safmin = dlamch( 'Safe minimum' )
290  safe1 = nz*safmin
291  safe2 = safe1 / eps
292 *
293 * Do for each right hand side
294 *
295  DO 250 j = 1, nrhs
296 *
297 * Compute residual R = B - op(A) * X,
298 * where op(A) = A, A**T, or A**H, depending on TRANS.
299 *
300  CALL zcopy( n, x( 1, j ), 1, work, 1 )
301  CALL ztrmv( uplo, trans, diag, n, a, lda, work, 1 )
302  CALL zaxpy( n, -one, b( 1, j ), 1, work, 1 )
303 *
304 * Compute componentwise relative backward error from formula
305 *
306 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
307 *
308 * where abs(Z) is the componentwise absolute value of the matrix
309 * or vector Z. If the i-th component of the denominator is less
310 * than SAFE2, then SAFE1 is added to the i-th components of the
311 * numerator and denominator before dividing.
312 *
313  DO 20 i = 1, n
314  rwork( i ) = cabs1( b( i, j ) )
315  20 CONTINUE
316 *
317  IF( notran ) THEN
318 *
319 * Compute abs(A)*abs(X) + abs(B).
320 *
321  IF( upper ) THEN
322  IF( nounit ) THEN
323  DO 40 k = 1, n
324  xk = cabs1( x( k, j ) )
325  DO 30 i = 1, k
326  rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
327  30 CONTINUE
328  40 CONTINUE
329  ELSE
330  DO 60 k = 1, n
331  xk = cabs1( x( k, j ) )
332  DO 50 i = 1, k - 1
333  rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
334  50 CONTINUE
335  rwork( k ) = rwork( k ) + xk
336  60 CONTINUE
337  END IF
338  ELSE
339  IF( nounit ) THEN
340  DO 80 k = 1, n
341  xk = cabs1( x( k, j ) )
342  DO 70 i = k, n
343  rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
344  70 CONTINUE
345  80 CONTINUE
346  ELSE
347  DO 100 k = 1, n
348  xk = cabs1( x( k, j ) )
349  DO 90 i = k + 1, n
350  rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
351  90 CONTINUE
352  rwork( k ) = rwork( k ) + xk
353  100 CONTINUE
354  END IF
355  END IF
356  ELSE
357 *
358 * Compute abs(A**H)*abs(X) + abs(B).
359 *
360  IF( upper ) THEN
361  IF( nounit ) THEN
362  DO 120 k = 1, n
363  s = zero
364  DO 110 i = 1, k
365  s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
366  110 CONTINUE
367  rwork( k ) = rwork( k ) + s
368  120 CONTINUE
369  ELSE
370  DO 140 k = 1, n
371  s = cabs1( x( k, j ) )
372  DO 130 i = 1, k - 1
373  s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
374  130 CONTINUE
375  rwork( k ) = rwork( k ) + s
376  140 CONTINUE
377  END IF
378  ELSE
379  IF( nounit ) THEN
380  DO 160 k = 1, n
381  s = zero
382  DO 150 i = k, n
383  s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
384  150 CONTINUE
385  rwork( k ) = rwork( k ) + s
386  160 CONTINUE
387  ELSE
388  DO 180 k = 1, n
389  s = cabs1( x( k, j ) )
390  DO 170 i = k + 1, n
391  s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
392  170 CONTINUE
393  rwork( k ) = rwork( k ) + s
394  180 CONTINUE
395  END IF
396  END IF
397  END IF
398  s = zero
399  DO 190 i = 1, n
400  IF( rwork( i ).GT.safe2 ) THEN
401  s = max( s, cabs1( work( i ) ) / rwork( i ) )
402  ELSE
403  s = max( s, ( cabs1( work( i ) )+safe1 ) /
404  $ ( rwork( i )+safe1 ) )
405  END IF
406  190 CONTINUE
407  berr( j ) = s
408 *
409 * Bound error from formula
410 *
411 * norm(X - XTRUE) / norm(X) .le. FERR =
412 * norm( abs(inv(op(A)))*
413 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
414 *
415 * where
416 * norm(Z) is the magnitude of the largest component of Z
417 * inv(op(A)) is the inverse of op(A)
418 * abs(Z) is the componentwise absolute value of the matrix or
419 * vector Z
420 * NZ is the maximum number of nonzeros in any row of A, plus 1
421 * EPS is machine epsilon
422 *
423 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
424 * is incremented by SAFE1 if the i-th component of
425 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
426 *
427 * Use ZLACN2 to estimate the infinity-norm of the matrix
428 * inv(op(A)) * diag(W),
429 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
430 *
431  DO 200 i = 1, n
432  IF( rwork( i ).GT.safe2 ) THEN
433  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
434  ELSE
435  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
436  $ safe1
437  END IF
438  200 CONTINUE
439 *
440  kase = 0
441  210 CONTINUE
442  CALL zlacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
443  IF( kase.NE.0 ) THEN
444  IF( kase.EQ.1 ) THEN
445 *
446 * Multiply by diag(W)*inv(op(A)**H).
447 *
448  CALL ztrsv( uplo, transt, diag, n, a, lda, work, 1 )
449  DO 220 i = 1, n
450  work( i ) = rwork( i )*work( i )
451  220 CONTINUE
452  ELSE
453 *
454 * Multiply by inv(op(A))*diag(W).
455 *
456  DO 230 i = 1, n
457  work( i ) = rwork( i )*work( i )
458  230 CONTINUE
459  CALL ztrsv( uplo, transn, diag, n, a, lda, work, 1 )
460  END IF
461  GO TO 210
462  END IF
463 *
464 * Normalize error.
465 *
466  lstres = zero
467  DO 240 i = 1, n
468  lstres = max( lstres, cabs1( x( i, j ) ) )
469  240 CONTINUE
470  IF( lstres.NE.zero )
471  $ ferr( j ) = ferr( j ) / lstres
472 *
473  250 CONTINUE
474 *
475  RETURN
476 *
477 * End of ZTRRFS
478 *
479  END
ztrsv
subroutine ztrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
ZTRSV
Definition: ztrsv.f:151
zcopy
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:135
ztrmv
subroutine ztrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
ZTRMV
Definition: ztrmv.f:149
ztrrfs
subroutine ztrrfs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
ZTRRFS
Definition: ztrrfs.f:184
zaxpy
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:53