LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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 complex16OTHERcomputational
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
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 ztpsv(UPLO, TRANS, DIAG, N, AP, X, INCX)
ZTPSV
Definition: ztpsv.f:144
subroutine ztpmv(UPLO, TRANS, DIAG, N, AP, X, INCX)
ZTPMV
Definition: ztpmv.f:142
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 ztprfs(UPLO, TRANS, DIAG, N, NRHS, AP, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
ZTPRFS
Definition: ztprfs.f:174