LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine cgtrfs ( character  TRANS,
integer  N,
integer  NRHS,
complex, dimension( * )  DL,
complex, dimension( * )  D,
complex, dimension( * )  DU,
complex, dimension( * )  DLF,
complex, dimension( * )  DF,
complex, dimension( * )  DUF,
complex, dimension( * )  DU2,
integer, dimension( * )  IPIV,
complex, dimension( ldb, * )  B,
integer  LDB,
complex, dimension( ldx, * )  X,
integer  LDX,
real, dimension( * )  FERR,
real, dimension( * )  BERR,
complex, dimension( * )  WORK,
real, dimension( * )  RWORK,
integer  INFO 
)

CGTRFS

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

Purpose: 
 CGTRFS 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 array, dimension (N-1)
          The (n-1) subdiagonal elements of A.
[in]D
          D is COMPLEX array, dimension (N)
          The diagonal elements of A.
[in]DU
          DU is COMPLEX array, dimension (N-1)
          The (n-1) superdiagonal elements of A.
[in]DLF
          DLF is COMPLEX array, dimension (N-1)
          The (n-1) multipliers that define the matrix L from the
          LU factorization of A as computed by CGTTRF.
[in]DF
          DF is COMPLEX array, dimension (N)
          The n diagonal elements of the upper triangular matrix U from
          the LU factorization of A.
[in]DUF
          DUF is COMPLEX array, dimension (N-1)
          The (n-1) elements of the first superdiagonal of U.
[in]DU2
          DU2 is COMPLEX 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 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 array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by CGTTRS.
          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 REAL 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 REAL 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 array, dimension (2*N)
[out]RWORK
          RWORK is REAL 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.
Date
September 2012

Definition at line 212 of file cgtrfs.f.

212 *
213 * -- LAPACK computational routine (version 3.4.2) --
214 * -- LAPACK is a software package provided by Univ. of Tennessee, --
215 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
216 * September 2012
217 *
218 * .. Scalar Arguments ..
219  CHARACTER trans
220  INTEGER info, ldb, ldx, n, nrhs
221 * ..
222 * .. Array Arguments ..
223  INTEGER ipiv( * )
224  REAL berr( * ), ferr( * ), rwork( * )
225  COMPLEX b( ldb, * ), d( * ), df( * ), dl( * ),
226  $ dlf( * ), du( * ), du2( * ), duf( * ),
227  $ work( * ), x( ldx, * )
228 * ..
229 *
230 * =====================================================================
231 *
232 * .. Parameters ..
233  INTEGER itmax
234  parameter ( itmax = 5 )
235  REAL zero, one
236  parameter ( zero = 0.0e+0, one = 1.0e+0 )
237  REAL two
238  parameter ( two = 2.0e+0 )
239  REAL three
240  parameter ( three = 3.0e+0 )
241 * ..
242 * .. Local Scalars ..
243  LOGICAL notran
244  CHARACTER transn, transt
245  INTEGER count, i, j, kase, nz
246  REAL eps, lstres, s, safe1, safe2, safmin
247  COMPLEX zdum
248 * ..
249 * .. Local Arrays ..
250  INTEGER isave( 3 )
251 * ..
252 * .. External Subroutines ..
253  EXTERNAL caxpy, ccopy, cgttrs, clacn2, clagtm, xerbla
254 * ..
255 * .. Intrinsic Functions ..
256  INTRINSIC abs, aimag, cmplx, max, real
257 * ..
258 * .. External Functions ..
259  LOGICAL lsame
260  REAL slamch
261  EXTERNAL lsame, slamch
262 * ..
263 * .. Statement Functions ..
264  REAL cabs1
265 * ..
266 * .. Statement Function definitions ..
267  cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( aimag( zdum ) )
268 * ..
269 * .. Executable Statements ..
270 *
271 * Test the input parameters.
272 *
273  info = 0
274  notran = lsame( trans, 'N' )
275  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
276  $ lsame( trans, 'C' ) ) THEN
277  info = -1
278  ELSE IF( n.LT.0 ) THEN
279  info = -2
280  ELSE IF( nrhs.LT.0 ) THEN
281  info = -3
282  ELSE IF( ldb.LT.max( 1, n ) ) THEN
283  info = -13
284  ELSE IF( ldx.LT.max( 1, n ) ) THEN
285  info = -15
286  END IF
287  IF( info.NE.0 ) THEN
288  CALL xerbla( 'CGTRFS', -info )
289  RETURN
290  END IF
291 *
292 * Quick return if possible
293 *
294  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
295  DO 10 j = 1, nrhs
296  ferr( j ) = zero
297  berr( j ) = zero
298  10 CONTINUE
299  RETURN
300  END IF
301 *
302  IF( notran ) THEN
303  transn = 'N'
304  transt = 'C'
305  ELSE
306  transn = 'C'
307  transt = 'N'
308  END IF
309 *
310 * NZ = maximum number of nonzero elements in each row of A, plus 1
311 *
312  nz = 4
313  eps = slamch( 'Epsilon' )
314  safmin = slamch( 'Safe minimum' )
315  safe1 = nz*safmin
316  safe2 = safe1 / eps
317 *
318 * Do for each right hand side
319 *
320  DO 110 j = 1, nrhs
321 *
322  count = 1
323  lstres = three
324  20 CONTINUE
325 *
326 * Loop until stopping criterion is satisfied.
327 *
328 * Compute residual R = B - op(A) * X,
329 * where op(A) = A, A**T, or A**H, depending on TRANS.
330 *
331  CALL ccopy( n, b( 1, j ), 1, work, 1 )
332  CALL clagtm( trans, n, 1, -one, dl, d, du, x( 1, j ), ldx, one,
333  $ work, n )
334 *
335 * Compute abs(op(A))*abs(x) + abs(b) for use in the backward
336 * error bound.
337 *
338  IF( notran ) THEN
339  IF( n.EQ.1 ) THEN
340  rwork( 1 ) = cabs1( b( 1, j ) ) +
341  $ cabs1( d( 1 ) )*cabs1( x( 1, j ) )
342  ELSE
343  rwork( 1 ) = cabs1( b( 1, j ) ) +
344  $ cabs1( d( 1 ) )*cabs1( x( 1, j ) ) +
345  $ cabs1( du( 1 ) )*cabs1( x( 2, j ) )
346  DO 30 i = 2, n - 1
347  rwork( i ) = cabs1( b( i, j ) ) +
348  $ cabs1( dl( i-1 ) )*cabs1( x( i-1, j ) ) +
349  $ cabs1( d( i ) )*cabs1( x( i, j ) ) +
350  $ cabs1( du( i ) )*cabs1( x( i+1, j ) )
351  30 CONTINUE
352  rwork( n ) = cabs1( b( n, j ) ) +
353  $ cabs1( dl( n-1 ) )*cabs1( x( n-1, j ) ) +
354  $ cabs1( d( n ) )*cabs1( x( n, j ) )
355  END IF
356  ELSE
357  IF( n.EQ.1 ) THEN
358  rwork( 1 ) = cabs1( b( 1, j ) ) +
359  $ cabs1( d( 1 ) )*cabs1( x( 1, j ) )
360  ELSE
361  rwork( 1 ) = cabs1( b( 1, j ) ) +
362  $ cabs1( d( 1 ) )*cabs1( x( 1, j ) ) +
363  $ cabs1( dl( 1 ) )*cabs1( x( 2, j ) )
364  DO 40 i = 2, n - 1
365  rwork( i ) = cabs1( b( i, j ) ) +
366  $ cabs1( du( i-1 ) )*cabs1( x( i-1, j ) ) +
367  $ cabs1( d( i ) )*cabs1( x( i, j ) ) +
368  $ cabs1( dl( i ) )*cabs1( x( i+1, j ) )
369  40 CONTINUE
370  rwork( n ) = cabs1( b( n, j ) ) +
371  $ cabs1( du( n-1 ) )*cabs1( x( n-1, j ) ) +
372  $ cabs1( d( n ) )*cabs1( x( n, j ) )
373  END IF
374  END IF
375 *
376 * Compute componentwise relative backward error from formula
377 *
378 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
379 *
380 * where abs(Z) is the componentwise absolute value of the matrix
381 * or vector Z. If the i-th component of the denominator is less
382 * than SAFE2, then SAFE1 is added to the i-th components of the
383 * numerator and denominator before dividing.
384 *
385  s = zero
386  DO 50 i = 1, n
387  IF( rwork( i ).GT.safe2 ) THEN
388  s = max( s, cabs1( work( i ) ) / rwork( i ) )
389  ELSE
390  s = max( s, ( cabs1( work( i ) )+safe1 ) /
391  $ ( rwork( i )+safe1 ) )
392  END IF
393  50 CONTINUE
394  berr( j ) = s
395 *
396 * Test stopping criterion. Continue iterating if
397 * 1) The residual BERR(J) is larger than machine epsilon, and
398 * 2) BERR(J) decreased by at least a factor of 2 during the
399 * last iteration, and
400 * 3) At most ITMAX iterations tried.
401 *
402  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
403  $ count.LE.itmax ) THEN
404 *
405 * Update solution and try again.
406 *
407  CALL cgttrs( trans, n, 1, dlf, df, duf, du2, ipiv, work, n,
408  $ info )
409  CALL caxpy( n, cmplx( one ), work, 1, x( 1, j ), 1 )
410  lstres = berr( j )
411  count = count + 1
412  GO TO 20
413  END IF
414 *
415 * Bound error from formula
416 *
417 * norm(X - XTRUE) / norm(X) .le. FERR =
418 * norm( abs(inv(op(A)))*
419 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
420 *
421 * where
422 * norm(Z) is the magnitude of the largest component of Z
423 * inv(op(A)) is the inverse of op(A)
424 * abs(Z) is the componentwise absolute value of the matrix or
425 * vector Z
426 * NZ is the maximum number of nonzeros in any row of A, plus 1
427 * EPS is machine epsilon
428 *
429 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
430 * is incremented by SAFE1 if the i-th component of
431 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
432 *
433 * Use CLACN2 to estimate the infinity-norm of the matrix
434 * inv(op(A)) * diag(W),
435 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
436 *
437  DO 60 i = 1, n
438  IF( rwork( i ).GT.safe2 ) THEN
439  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
440  ELSE
441  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
442  $ safe1
443  END IF
444  60 CONTINUE
445 *
446  kase = 0
447  70 CONTINUE
448  CALL clacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
449  IF( kase.NE.0 ) THEN
450  IF( kase.EQ.1 ) THEN
451 *
452 * Multiply by diag(W)*inv(op(A)**H).
453 *
454  CALL cgttrs( transt, n, 1, dlf, df, duf, du2, ipiv, 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 cgttrs( transn, n, 1, dlf, df, duf, du2, ipiv, work,
467  $ n, info )
468  END IF
469  GO TO 70
470  END IF
471 *
472 * Normalize error.
473 *
474  lstres = zero
475  DO 100 i = 1, n
476  lstres = max( lstres, cabs1( x( i, j ) ) )
477  100 CONTINUE
478  IF( lstres.NE.zero )
479  $ ferr( j ) = ferr( j ) / lstres
480 *
481  110 CONTINUE
482 *
483  RETURN
484 *
485 * End of CGTRFS
486 *
cgttrs
subroutine cgttrs(TRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB, INFO)
CGTTRS
Definition: cgttrs.f:140
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
clagtm
subroutine clagtm(TRANS, N, NRHS, ALPHA, DL, D, DU, X, LDX, BETA, B, LDB)
CLAGTM performs a matrix-matrix product of the form C = αAB+βC, where A is a tridiagonal matrix...
Definition: clagtm.f:147
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
caxpy
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:53
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
clacn2
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:135

Here is the call graph for this function:

Here is the caller graph for this function: