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

◆ 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
cgttrs
subroutine cgttrs(trans, n, nrhs, dl, d, du, du2, ipiv, b, ldb, info)
CGTTRS
Definition cgttrs.f:138
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:133
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:145
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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: