LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cgbrfs.f
Go to the documentation of this file.
1*> \brief \b CGBRFS
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CGBRFS + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgbrfs.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgbrfs.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgbrfs.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CGBRFS( 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* REAL BERR( * ), FERR( * ), RWORK( * )
32* COMPLEX AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
33* $ WORK( * ), X( LDX, * )
34* ..
35*
36*
37*> \par Purpose:
38* =============
39*>
40*> \verbatim
41*>
42*> CGBRFS 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 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 array, dimension (LDAFB,N)
102*> Details of the LU factorization of the band matrix A, as
103*> computed by CGBTRF. 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 CGBTRF; 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 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 array, dimension (LDX,NRHS)
137*> On entry, the solution matrix X, as computed by CGBTRS.
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 REAL 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 REAL 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 array, dimension (2*N)
171*> \endverbatim
172*>
173*> \param[out] RWORK
174*> \verbatim
175*> RWORK is REAL 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 gbrfs
201*
202* =====================================================================
203 SUBROUTINE cgbrfs( 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 REAL BERR( * ), FERR( * ), RWORK( * )
218 COMPLEX AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
219 $ work( * ), x( ldx, * )
220* ..
221*
222* =====================================================================
223*
224* .. Parameters ..
225 INTEGER ITMAX
226 PARAMETER ( ITMAX = 5 )
227 REAL ZERO
228 parameter( zero = 0.0e+0 )
229 COMPLEX CONE
230 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
231 REAL TWO
232 parameter( two = 2.0e+0 )
233 REAL THREE
234 parameter( three = 3.0e+0 )
235* ..
236* .. Local Scalars ..
237 LOGICAL NOTRAN
238 CHARACTER TRANSN, TRANST
239 INTEGER COUNT, I, J, K, KASE, KK, NZ
240 REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
241 COMPLEX ZDUM
242* ..
243* .. Local Arrays ..
244 INTEGER ISAVE( 3 )
245* ..
246* .. External Subroutines ..
247 EXTERNAL caxpy, ccopy, cgbmv, cgbtrs, clacn2, xerbla
248* ..
249* .. Intrinsic Functions ..
250 INTRINSIC abs, aimag, max, min, real
251* ..
252* .. External Functions ..
253 LOGICAL LSAME
254 REAL SLAMCH
255 EXTERNAL lsame, slamch
256* ..
257* .. Statement Functions ..
258 REAL CABS1
259* ..
260* .. Statement Function definitions ..
261 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( 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( 'CGBRFS', -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 = slamch( 'Epsilon' )
316 safmin = slamch( '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 ccopy( n, b( 1, j ), 1, work, 1 )
334 CALL cgbmv( 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 cgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv, work, n,
393 $ info )
394 CALL caxpy( 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 CLACN2 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 clacn2( 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 cgbtrs( 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 cgbtrs( 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 CGBRFS
471*
472 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgbmv(trans, m, n, kl, ku, alpha, a, lda, x, incx, beta, y, incy)
CGBMV
Definition cgbmv.f:190
subroutine cgbrfs(trans, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv, b, ldb, x, ldx, ferr, berr, work, rwork, info)
CGBRFS
Definition cgbrfs.f:206
subroutine cgbtrs(trans, n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info)
CGBTRS
Definition cgbtrs.f:138
subroutine clacn2(n, v, x, est, kase, isave)
CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition clacn2.f:133