LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cgtrfs()

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.

Definition at line 207 of file cgtrfs.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  REAL BERR( * ), FERR( * ), RWORK( * )
222  COMPLEX 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  REAL ZERO, ONE
233  parameter( zero = 0.0e+0, one = 1.0e+0 )
234  REAL TWO
235  parameter( two = 2.0e+0 )
236  REAL THREE
237  parameter( three = 3.0e+0 )
238 * ..
239 * .. Local Scalars ..
240  LOGICAL NOTRAN
241  CHARACTER TRANSN, TRANST
242  INTEGER COUNT, I, J, KASE, NZ
243  REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
244  COMPLEX ZDUM
245 * ..
246 * .. Local Arrays ..
247  INTEGER ISAVE( 3 )
248 * ..
249 * .. External Subroutines ..
250  EXTERNAL caxpy, ccopy, cgttrs, clacn2, clagtm, xerbla
251 * ..
252 * .. Intrinsic Functions ..
253  INTRINSIC abs, aimag, cmplx, max, real
254 * ..
255 * .. External Functions ..
256  LOGICAL LSAME
257  REAL SLAMCH
258  EXTERNAL lsame, slamch
259 * ..
260 * .. Statement Functions ..
261  REAL CABS1
262 * ..
263 * .. Statement Function definitions ..
264  cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( 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( 'CGTRFS', -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 = slamch( 'Epsilon' )
311  safmin = slamch( '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 ccopy( n, b( 1, j ), 1, work, 1 )
329  CALL clagtm( 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 cgttrs( trans, n, 1, dlf, df, duf, du2, ipiv, work, n,
405  $ info )
406  CALL caxpy( n, cmplx( 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 CLACN2 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 clacn2( 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 cgttrs( 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 cgttrs( 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 CGTRFS
483 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
subroutine cgttrs(TRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB, INFO)
CGTTRS
Definition: cgttrs.f:138
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:145
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
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: