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