LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zgtrfs()

subroutine zgtrfs ( character  trans,
integer  n,
integer  nrhs,
complex*16, dimension( * )  dl,
complex*16, dimension( * )  d,
complex*16, dimension( * )  du,
complex*16, dimension( * )  dlf,
complex*16, dimension( * )  df,
complex*16, dimension( * )  duf,
complex*16, dimension( * )  du2,
integer, dimension( * )  ipiv,
complex*16, dimension( ldb, * )  b,
integer  ldb,
complex*16, dimension( ldx, * )  x,
integer  ldx,
double precision, dimension( * )  ferr,
double precision, dimension( * )  berr,
complex*16, dimension( * )  work,
double precision, dimension( * )  rwork,
integer  info 
)

ZGTRFS

Download ZGTRFS + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGTRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is tridiagonal, and provides
 error bounds and backward error estimates for the solution.
Parameters
[in]TRANS
          TRANS is CHARACTER*1
          Specifies the form of the system of equations:
          = 'N':  A * X = B     (No transpose)
          = 'T':  A**T * X = B  (Transpose)
          = 'C':  A**H * X = B  (Conjugate transpose)
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrix B.  NRHS >= 0.
[in]DL
          DL is COMPLEX*16 array, dimension (N-1)
          The (n-1) subdiagonal elements of A.
[in]D
          D is COMPLEX*16 array, dimension (N)
          The diagonal elements of A.
[in]DU
          DU is COMPLEX*16 array, dimension (N-1)
          The (n-1) superdiagonal elements of A.
[in]DLF
          DLF is COMPLEX*16 array, dimension (N-1)
          The (n-1) multipliers that define the matrix L from the
          LU factorization of A as computed by ZGTTRF.
[in]DF
          DF is COMPLEX*16 array, dimension (N)
          The n diagonal elements of the upper triangular matrix U from
          the LU factorization of A.
[in]DUF
          DUF is COMPLEX*16 array, dimension (N-1)
          The (n-1) elements of the first superdiagonal of U.
[in]DU2
          DU2 is COMPLEX*16 array, dimension (N-2)
          The (n-2) elements of the second superdiagonal of U.
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          The pivot indices; for 1 <= i <= n, row i of the matrix was
          interchanged with row IPIV(i).  IPIV(i) will always be either
          i or i+1; IPIV(i) = i indicates a row interchange was not
          required.
[in]B
          B is COMPLEX*16 array, dimension (LDB,NRHS)
          The right hand side matrix B.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[in,out]X
          X is COMPLEX*16 array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by ZGTTRS.
          On exit, the improved solution matrix X.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.  LDX >= max(1,N).
[out]FERR
          FERR is DOUBLE PRECISION array, dimension (NRHS)
          The estimated forward error bound for each solution vector
          X(j) (the j-th column of the solution matrix X).
          If XTRUE is the true solution corresponding to X(j), FERR(j)
          is an estimated upper bound for the magnitude of the largest
          element in (X(j) - XTRUE) divided by the magnitude of the
          largest element in X(j).  The estimate is as reliable as
          the estimate for RCOND, and is almost always a slight
          overestimate of the true error.
[out]BERR
          BERR is DOUBLE PRECISION array, dimension (NRHS)
          The componentwise relative backward error of each solution
          vector X(j) (i.e., the smallest relative change in
          any element of A or B that makes X(j) an exact solution).
[out]WORK
          WORK is COMPLEX*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Internal Parameters:
  ITMAX is the maximum number of steps of iterative refinement.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 207 of file zgtrfs.f.

210*
211* -- LAPACK computational routine --
212* -- LAPACK is a software package provided by Univ. of Tennessee, --
213* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
214*
215* .. Scalar Arguments ..
216 CHARACTER TRANS
217 INTEGER INFO, LDB, LDX, N, NRHS
218* ..
219* .. Array Arguments ..
220 INTEGER IPIV( * )
221 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
222 COMPLEX*16 B( LDB, * ), D( * ), DF( * ), DL( * ),
223 $ DLF( * ), DU( * ), DU2( * ), DUF( * ),
224 $ WORK( * ), X( LDX, * )
225* ..
226*
227* =====================================================================
228*
229* .. Parameters ..
230 INTEGER ITMAX
231 parameter( itmax = 5 )
232 DOUBLE PRECISION ZERO, ONE
233 parameter( zero = 0.0d+0, one = 1.0d+0 )
234 DOUBLE PRECISION TWO
235 parameter( two = 2.0d+0 )
236 DOUBLE PRECISION THREE
237 parameter( three = 3.0d+0 )
238* ..
239* .. Local Scalars ..
240 LOGICAL NOTRAN
241 CHARACTER TRANSN, TRANST
242 INTEGER COUNT, I, J, KASE, NZ
243 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
244 COMPLEX*16 ZDUM
245* ..
246* .. Local Arrays ..
247 INTEGER ISAVE( 3 )
248* ..
249* .. External Subroutines ..
250 EXTERNAL xerbla, zaxpy, zcopy, zgttrs, zlacn2, 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, one,
330 $ work, n )
331*
332* Compute abs(op(A))*abs(x) + abs(b) for use in the backward
333* error bound.
334*
335 IF( notran ) THEN
336 IF( n.EQ.1 ) THEN
337 rwork( 1 ) = cabs1( b( 1, j ) ) +
338 $ cabs1( d( 1 ) )*cabs1( x( 1, j ) )
339 ELSE
340 rwork( 1 ) = cabs1( b( 1, j ) ) +
341 $ cabs1( d( 1 ) )*cabs1( x( 1, j ) ) +
342 $ cabs1( du( 1 ) )*cabs1( x( 2, j ) )
343 DO 30 i = 2, n - 1
344 rwork( i ) = cabs1( b( i, j ) ) +
345 $ cabs1( dl( i-1 ) )*cabs1( x( i-1, j ) ) +
346 $ cabs1( d( i ) )*cabs1( x( i, j ) ) +
347 $ cabs1( du( i ) )*cabs1( x( i+1, j ) )
348 30 CONTINUE
349 rwork( n ) = cabs1( b( n, j ) ) +
350 $ cabs1( dl( n-1 ) )*cabs1( x( n-1, j ) ) +
351 $ cabs1( d( n ) )*cabs1( x( n, j ) )
352 END IF
353 ELSE
354 IF( n.EQ.1 ) THEN
355 rwork( 1 ) = cabs1( b( 1, j ) ) +
356 $ cabs1( d( 1 ) )*cabs1( x( 1, j ) )
357 ELSE
358 rwork( 1 ) = cabs1( b( 1, j ) ) +
359 $ cabs1( d( 1 ) )*cabs1( x( 1, j ) ) +
360 $ cabs1( dl( 1 ) )*cabs1( x( 2, j ) )
361 DO 40 i = 2, n - 1
362 rwork( i ) = cabs1( b( i, j ) ) +
363 $ cabs1( du( i-1 ) )*cabs1( x( i-1, j ) ) +
364 $ cabs1( d( i ) )*cabs1( x( i, j ) ) +
365 $ cabs1( dl( i ) )*cabs1( x( i+1, j ) )
366 40 CONTINUE
367 rwork( n ) = cabs1( b( n, j ) ) +
368 $ cabs1( du( n-1 ) )*cabs1( x( n-1, j ) ) +
369 $ cabs1( d( n ) )*cabs1( x( n, j ) )
370 END IF
371 END IF
372*
373* Compute componentwise relative backward error from formula
374*
375* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
376*
377* where abs(Z) is the componentwise absolute value of the matrix
378* or vector Z. If the i-th component of the denominator is less
379* than SAFE2, then SAFE1 is added to the i-th components of the
380* numerator and denominator before dividing.
381*
382 s = zero
383 DO 50 i = 1, n
384 IF( rwork( i ).GT.safe2 ) THEN
385 s = max( s, cabs1( work( i ) ) / rwork( i ) )
386 ELSE
387 s = max( s, ( cabs1( work( i ) )+safe1 ) /
388 $ ( rwork( i )+safe1 ) )
389 END IF
390 50 CONTINUE
391 berr( j ) = s
392*
393* Test stopping criterion. Continue iterating if
394* 1) The residual BERR(J) is larger than machine epsilon, and
395* 2) BERR(J) decreased by at least a factor of 2 during the
396* last iteration, and
397* 3) At most ITMAX iterations tried.
398*
399 IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
400 $ count.LE.itmax ) THEN
401*
402* Update solution and try again.
403*
404 CALL zgttrs( trans, n, 1, dlf, df, duf, du2, ipiv, work, n,
405 $ info )
406 CALL zaxpy( n, dcmplx( one ), work, 1, x( 1, j ), 1 )
407 lstres = berr( j )
408 count = count + 1
409 GO TO 20
410 END IF
411*
412* Bound error from formula
413*
414* norm(X - XTRUE) / norm(X) .le. FERR =
415* norm( abs(inv(op(A)))*
416* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
417*
418* where
419* norm(Z) is the magnitude of the largest component of Z
420* inv(op(A)) is the inverse of op(A)
421* abs(Z) is the componentwise absolute value of the matrix or
422* vector Z
423* NZ is the maximum number of nonzeros in any row of A, plus 1
424* EPS is machine epsilon
425*
426* The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
427* is incremented by SAFE1 if the i-th component of
428* abs(op(A))*abs(X) + abs(B) is less than SAFE2.
429*
430* Use ZLACN2 to estimate the infinity-norm of the matrix
431* inv(op(A)) * diag(W),
432* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
433*
434 DO 60 i = 1, n
435 IF( rwork( i ).GT.safe2 ) THEN
436 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
437 ELSE
438 rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
439 $ safe1
440 END IF
441 60 CONTINUE
442*
443 kase = 0
444 70 CONTINUE
445 CALL zlacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
446 IF( kase.NE.0 ) THEN
447 IF( kase.EQ.1 ) THEN
448*
449* Multiply by diag(W)*inv(op(A)**H).
450*
451 CALL zgttrs( transt, n, 1, dlf, df, duf, du2, ipiv, work,
452 $ n, info )
453 DO 80 i = 1, n
454 work( i ) = rwork( i )*work( i )
455 80 CONTINUE
456 ELSE
457*
458* Multiply by inv(op(A))*diag(W).
459*
460 DO 90 i = 1, n
461 work( i ) = rwork( i )*work( i )
462 90 CONTINUE
463 CALL zgttrs( transn, n, 1, dlf, df, duf, du2, ipiv, work,
464 $ n, info )
465 END IF
466 GO TO 70
467 END IF
468*
469* Normalize error.
470*
471 lstres = zero
472 DO 100 i = 1, n
473 lstres = max( lstres, cabs1( x( i, j ) ) )
474 100 CONTINUE
475 IF( lstres.NE.zero )
476 $ ferr( j ) = ferr( j ) / lstres
477*
478 110 CONTINUE
479*
480 RETURN
481*
482* End of ZGTRFS
483*
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
zgttrs
subroutine zgttrs(trans, n, nrhs, dl, d, du, du2, ipiv, b, ldb, info)
ZGTTRS
Definition zgttrs.f:138
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:133
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:145
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the call graph for this function:
Here is the caller graph for this function: