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