LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
zgbrfs.f
Go to the documentation of this file.
1 *> \brief \b ZGBRFS
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZGBRFS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgbrfs.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgbrfs.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgbrfs.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZGBRFS( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB,
22 * IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK,
23 * INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER TRANS
27 * INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS
28 * ..
29 * .. Array Arguments ..
30 * INTEGER IPIV( * )
31 * DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
32 * COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
33 * $ WORK( * ), X( LDX, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> ZGBRFS improves the computed solution to a system of linear
43 *> equations when the coefficient matrix is banded, and provides
44 *> error bounds and backward error estimates for the solution.
45 *> \endverbatim
46 *
47 * Arguments:
48 * ==========
49 *
50 *> \param[in] TRANS
51 *> \verbatim
52 *> TRANS is CHARACTER*1
53 *> Specifies the form of the system of equations:
54 *> = 'N': A * X = B (No transpose)
55 *> = 'T': A**T * X = B (Transpose)
56 *> = 'C': A**H * X = B (Conjugate transpose)
57 *> \endverbatim
58 *>
59 *> \param[in] N
60 *> \verbatim
61 *> N is INTEGER
62 *> The order of the matrix A. N >= 0.
63 *> \endverbatim
64 *>
65 *> \param[in] KL
66 *> \verbatim
67 *> KL is INTEGER
68 *> The number of subdiagonals within the band of A. KL >= 0.
69 *> \endverbatim
70 *>
71 *> \param[in] KU
72 *> \verbatim
73 *> KU is INTEGER
74 *> The number of superdiagonals within the band of A. KU >= 0.
75 *> \endverbatim
76 *>
77 *> \param[in] NRHS
78 *> \verbatim
79 *> NRHS is INTEGER
80 *> The number of right hand sides, i.e., the number of columns
81 *> of the matrices B and X. NRHS >= 0.
82 *> \endverbatim
83 *>
84 *> \param[in] AB
85 *> \verbatim
86 *> AB is COMPLEX*16 array, dimension (LDAB,N)
87 *> The original band matrix A, stored in rows 1 to KL+KU+1.
88 *> The j-th column of A is stored in the j-th column of the
89 *> array AB as follows:
90 *> AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl).
91 *> \endverbatim
92 *>
93 *> \param[in] LDAB
94 *> \verbatim
95 *> LDAB is INTEGER
96 *> The leading dimension of the array AB. LDAB >= KL+KU+1.
97 *> \endverbatim
98 *>
99 *> \param[in] AFB
100 *> \verbatim
101 *> AFB is COMPLEX*16 array, dimension (LDAFB,N)
102 *> Details of the LU factorization of the band matrix A, as
103 *> computed by ZGBTRF. U is stored as an upper triangular band
104 *> matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
105 *> the multipliers used during the factorization are stored in
106 *> rows KL+KU+2 to 2*KL+KU+1.
107 *> \endverbatim
108 *>
109 *> \param[in] LDAFB
110 *> \verbatim
111 *> LDAFB is INTEGER
112 *> The leading dimension of the array AFB. LDAFB >= 2*KL*KU+1.
113 *> \endverbatim
114 *>
115 *> \param[in] IPIV
116 *> \verbatim
117 *> IPIV is INTEGER array, dimension (N)
118 *> The pivot indices from ZGBTRF; for 1<=i<=N, row i of the
119 *> matrix was interchanged with row IPIV(i).
120 *> \endverbatim
121 *>
122 *> \param[in] B
123 *> \verbatim
124 *> B is COMPLEX*16 array, dimension (LDB,NRHS)
125 *> The right hand side matrix B.
126 *> \endverbatim
127 *>
128 *> \param[in] LDB
129 *> \verbatim
130 *> LDB is INTEGER
131 *> The leading dimension of the array B. LDB >= max(1,N).
132 *> \endverbatim
133 *>
134 *> \param[in,out] X
135 *> \verbatim
136 *> X is COMPLEX*16 array, dimension (LDX,NRHS)
137 *> On entry, the solution matrix X, as computed by ZGBTRS.
138 *> On exit, the improved solution matrix X.
139 *> \endverbatim
140 *>
141 *> \param[in] LDX
142 *> \verbatim
143 *> LDX is INTEGER
144 *> The leading dimension of the array X. LDX >= max(1,N).
145 *> \endverbatim
146 *>
147 *> \param[out] FERR
148 *> \verbatim
149 *> FERR is DOUBLE PRECISION array, dimension (NRHS)
150 *> The estimated forward error bound for each solution vector
151 *> X(j) (the j-th column of the solution matrix X).
152 *> If XTRUE is the true solution corresponding to X(j), FERR(j)
153 *> is an estimated upper bound for the magnitude of the largest
154 *> element in (X(j) - XTRUE) divided by the magnitude of the
155 *> largest element in X(j). The estimate is as reliable as
156 *> the estimate for RCOND, and is almost always a slight
157 *> overestimate of the true error.
158 *> \endverbatim
159 *>
160 *> \param[out] BERR
161 *> \verbatim
162 *> BERR is DOUBLE PRECISION array, dimension (NRHS)
163 *> The componentwise relative backward error of each solution
164 *> vector X(j) (i.e., the smallest relative change in
165 *> any element of A or B that makes X(j) an exact solution).
166 *> \endverbatim
167 *>
168 *> \param[out] WORK
169 *> \verbatim
170 *> WORK is COMPLEX*16 array, dimension (2*N)
171 *> \endverbatim
172 *>
173 *> \param[out] RWORK
174 *> \verbatim
175 *> RWORK is DOUBLE PRECISION array, dimension (N)
176 *> \endverbatim
177 *>
178 *> \param[out] INFO
179 *> \verbatim
180 *> INFO is INTEGER
181 *> = 0: successful exit
182 *> < 0: if INFO = -i, the i-th argument had an illegal value
183 *> \endverbatim
184 *
185 *> \par Internal Parameters:
186 * =========================
187 *>
188 *> \verbatim
189 *> ITMAX is the maximum number of steps of iterative refinement.
190 *> \endverbatim
191 *
192 * Authors:
193 * ========
194 *
195 *> \author Univ. of Tennessee
196 *> \author Univ. of California Berkeley
197 *> \author Univ. of Colorado Denver
198 *> \author NAG Ltd.
199 *
200 *> \ingroup complex16GBcomputational
201 *
202 * =====================================================================
203  SUBROUTINE zgbrfs( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB,
204  $ IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK,
205  $ INFO )
206 *
207 * -- LAPACK computational routine --
208 * -- LAPACK is a software package provided by Univ. of Tennessee, --
209 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
210 *
211 * .. Scalar Arguments ..
212  CHARACTER TRANS
213  INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS
214 * ..
215 * .. Array Arguments ..
216  INTEGER IPIV( * )
217  DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
218  COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
219  $ work( * ), x( ldx, * )
220 * ..
221 *
222 * =====================================================================
223 *
224 * .. Parameters ..
225  INTEGER ITMAX
226  PARAMETER ( ITMAX = 5 )
227  DOUBLE PRECISION ZERO
228  parameter( zero = 0.0d+0 )
229  COMPLEX*16 CONE
230  parameter( cone = ( 1.0d+0, 0.0d+0 ) )
231  DOUBLE PRECISION TWO
232  parameter( two = 2.0d+0 )
233  DOUBLE PRECISION THREE
234  parameter( three = 3.0d+0 )
235 * ..
236 * .. Local Scalars ..
237  LOGICAL NOTRAN
238  CHARACTER TRANSN, TRANST
239  INTEGER COUNT, I, J, K, KASE, KK, NZ
240  DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
241  COMPLEX*16 ZDUM
242 * ..
243 * .. Local Arrays ..
244  INTEGER ISAVE( 3 )
245 * ..
246 * .. External Subroutines ..
247  EXTERNAL xerbla, zaxpy, zcopy, zgbmv, zgbtrs, zlacn2
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC abs, dble, dimag, max, min
251 * ..
252 * .. External Functions ..
253  LOGICAL LSAME
254  DOUBLE PRECISION DLAMCH
255  EXTERNAL lsame, dlamch
256 * ..
257 * .. Statement Functions ..
258  DOUBLE PRECISION CABS1
259 * ..
260 * .. Statement Function definitions ..
261  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
262 * ..
263 * .. Executable Statements ..
264 *
265 * Test the input parameters.
266 *
267  info = 0
268  notran = lsame( trans, 'N' )
269  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
270  $ lsame( trans, 'C' ) ) THEN
271  info = -1
272  ELSE IF( n.LT.0 ) THEN
273  info = -2
274  ELSE IF( kl.LT.0 ) THEN
275  info = -3
276  ELSE IF( ku.LT.0 ) THEN
277  info = -4
278  ELSE IF( nrhs.LT.0 ) THEN
279  info = -5
280  ELSE IF( ldab.LT.kl+ku+1 ) THEN
281  info = -7
282  ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
283  info = -9
284  ELSE IF( ldb.LT.max( 1, n ) ) THEN
285  info = -12
286  ELSE IF( ldx.LT.max( 1, n ) ) THEN
287  info = -14
288  END IF
289  IF( info.NE.0 ) THEN
290  CALL xerbla( 'ZGBRFS', -info )
291  RETURN
292  END IF
293 *
294 * Quick return if possible
295 *
296  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
297  DO 10 j = 1, nrhs
298  ferr( j ) = zero
299  berr( j ) = zero
300  10 CONTINUE
301  RETURN
302  END IF
303 *
304  IF( notran ) THEN
305  transn = 'N'
306  transt = 'C'
307  ELSE
308  transn = 'C'
309  transt = 'N'
310  END IF
311 *
312 * NZ = maximum number of nonzero elements in each row of A, plus 1
313 *
314  nz = min( kl+ku+2, n+1 )
315  eps = dlamch( 'Epsilon' )
316  safmin = dlamch( 'Safe minimum' )
317  safe1 = nz*safmin
318  safe2 = safe1 / eps
319 *
320 * Do for each right hand side
321 *
322  DO 140 j = 1, nrhs
323 *
324  count = 1
325  lstres = three
326  20 CONTINUE
327 *
328 * Loop until stopping criterion is satisfied.
329 *
330 * Compute residual R = B - op(A) * X,
331 * where op(A) = A, A**T, or A**H, depending on TRANS.
332 *
333  CALL zcopy( n, b( 1, j ), 1, work, 1 )
334  CALL zgbmv( trans, n, n, kl, ku, -cone, ab, ldab, x( 1, j ), 1,
335  $ cone, work, 1 )
336 *
337 * Compute componentwise relative backward error from formula
338 *
339 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
340 *
341 * where abs(Z) is the componentwise absolute value of the matrix
342 * or vector Z. If the i-th component of the denominator is less
343 * than SAFE2, then SAFE1 is added to the i-th components of the
344 * numerator and denominator before dividing.
345 *
346  DO 30 i = 1, n
347  rwork( i ) = cabs1( b( i, j ) )
348  30 CONTINUE
349 *
350 * Compute abs(op(A))*abs(X) + abs(B).
351 *
352  IF( notran ) THEN
353  DO 50 k = 1, n
354  kk = ku + 1 - k
355  xk = cabs1( x( k, j ) )
356  DO 40 i = max( 1, k-ku ), min( n, k+kl )
357  rwork( i ) = rwork( i ) + cabs1( ab( kk+i, k ) )*xk
358  40 CONTINUE
359  50 CONTINUE
360  ELSE
361  DO 70 k = 1, n
362  s = zero
363  kk = ku + 1 - k
364  DO 60 i = max( 1, k-ku ), min( n, k+kl )
365  s = s + cabs1( ab( kk+i, k ) )*cabs1( x( i, j ) )
366  60 CONTINUE
367  rwork( k ) = rwork( k ) + s
368  70 CONTINUE
369  END IF
370  s = zero
371  DO 80 i = 1, n
372  IF( rwork( i ).GT.safe2 ) THEN
373  s = max( s, cabs1( work( i ) ) / rwork( i ) )
374  ELSE
375  s = max( s, ( cabs1( work( i ) )+safe1 ) /
376  $ ( rwork( i )+safe1 ) )
377  END IF
378  80 CONTINUE
379  berr( j ) = s
380 *
381 * Test stopping criterion. Continue iterating if
382 * 1) The residual BERR(J) is larger than machine epsilon, and
383 * 2) BERR(J) decreased by at least a factor of 2 during the
384 * last iteration, and
385 * 3) At most ITMAX iterations tried.
386 *
387  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
388  $ count.LE.itmax ) THEN
389 *
390 * Update solution and try again.
391 *
392  CALL zgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv, work, n,
393  $ info )
394  CALL zaxpy( n, cone, work, 1, x( 1, j ), 1 )
395  lstres = berr( j )
396  count = count + 1
397  GO TO 20
398  END IF
399 *
400 * Bound error from formula
401 *
402 * norm(X - XTRUE) / norm(X) .le. FERR =
403 * norm( abs(inv(op(A)))*
404 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
405 *
406 * where
407 * norm(Z) is the magnitude of the largest component of Z
408 * inv(op(A)) is the inverse of op(A)
409 * abs(Z) is the componentwise absolute value of the matrix or
410 * vector Z
411 * NZ is the maximum number of nonzeros in any row of A, plus 1
412 * EPS is machine epsilon
413 *
414 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
415 * is incremented by SAFE1 if the i-th component of
416 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
417 *
418 * Use ZLACN2 to estimate the infinity-norm of the matrix
419 * inv(op(A)) * diag(W),
420 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
421 *
422  DO 90 i = 1, n
423  IF( rwork( i ).GT.safe2 ) THEN
424  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
425  ELSE
426  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
427  $ safe1
428  END IF
429  90 CONTINUE
430 *
431  kase = 0
432  100 CONTINUE
433  CALL zlacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
434  IF( kase.NE.0 ) THEN
435  IF( kase.EQ.1 ) THEN
436 *
437 * Multiply by diag(W)*inv(op(A)**H).
438 *
439  CALL zgbtrs( transt, n, kl, ku, 1, afb, ldafb, ipiv,
440  $ work, n, info )
441  DO 110 i = 1, n
442  work( i ) = rwork( i )*work( i )
443  110 CONTINUE
444  ELSE
445 *
446 * Multiply by inv(op(A))*diag(W).
447 *
448  DO 120 i = 1, n
449  work( i ) = rwork( i )*work( i )
450  120 CONTINUE
451  CALL zgbtrs( transn, n, kl, ku, 1, afb, ldafb, ipiv,
452  $ work, n, info )
453  END IF
454  GO TO 100
455  END IF
456 *
457 * Normalize error.
458 *
459  lstres = zero
460  DO 130 i = 1, n
461  lstres = max( lstres, cabs1( x( i, j ) ) )
462  130 CONTINUE
463  IF( lstres.NE.zero )
464  $ ferr( j ) = ferr( j ) / lstres
465 *
466  140 CONTINUE
467 *
468  RETURN
469 *
470 * End of ZGBRFS
471 *
472  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 zgbmv(TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGBMV
Definition: zgbmv.f:187
subroutine zgbrfs(TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
ZGBRFS
Definition: zgbrfs.f:206
subroutine zgbtrs(TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)
ZGBTRS
Definition: zgbtrs.f:138
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