LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ zptrfs()

subroutine zptrfs ( character  UPLO,
integer  N,
integer  NRHS,
double precision, dimension( * )  D,
complex*16, dimension( * )  E,
double precision, dimension( * )  DF,
complex*16, dimension( * )  EF,
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 
)

ZPTRFS

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

Purpose:
 ZPTRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is Hermitian positive definite
 and tridiagonal, and provides error bounds and backward error
 estimates for the solution.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the superdiagonal or the subdiagonal of the
          tridiagonal matrix A is stored and the form of the
          factorization:
          = 'U':  E is the superdiagonal of A, and A = U**H*D*U;
          = 'L':  E is the subdiagonal of A, and A = L*D*L**H.
          (The two forms are equivalent if A is real.)
[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]D
          D is DOUBLE PRECISION array, dimension (N)
          The n real diagonal elements of the tridiagonal matrix A.
[in]E
          E is COMPLEX*16 array, dimension (N-1)
          The (n-1) off-diagonal elements of the tridiagonal matrix A
          (see UPLO).
[in]DF
          DF is DOUBLE PRECISION array, dimension (N)
          The n diagonal elements of the diagonal matrix D from
          the factorization computed by ZPTTRF.
[in]EF
          EF is COMPLEX*16 array, dimension (N-1)
          The (n-1) off-diagonal elements of the unit bidiagonal
          factor U or L from the factorization computed by ZPTTRF
          (see UPLO).
[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 ZPTTRS.
          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 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).
[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 (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 181 of file zptrfs.f.

183 *
184 * -- LAPACK computational routine --
185 * -- LAPACK is a software package provided by Univ. of Tennessee, --
186 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
187 *
188 * .. Scalar Arguments ..
189  CHARACTER UPLO
190  INTEGER INFO, LDB, LDX, N, NRHS
191 * ..
192 * .. Array Arguments ..
193  DOUBLE PRECISION BERR( * ), D( * ), DF( * ), FERR( * ),
194  $ RWORK( * )
195  COMPLEX*16 B( LDB, * ), E( * ), EF( * ), WORK( * ),
196  $ X( LDX, * )
197 * ..
198 *
199 * =====================================================================
200 *
201 * .. Parameters ..
202  INTEGER ITMAX
203  parameter( itmax = 5 )
204  DOUBLE PRECISION ZERO
205  parameter( zero = 0.0d+0 )
206  DOUBLE PRECISION ONE
207  parameter( one = 1.0d+0 )
208  DOUBLE PRECISION TWO
209  parameter( two = 2.0d+0 )
210  DOUBLE PRECISION THREE
211  parameter( three = 3.0d+0 )
212 * ..
213 * .. Local Scalars ..
214  LOGICAL UPPER
215  INTEGER COUNT, I, IX, J, NZ
216  DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
217  COMPLEX*16 BI, CX, DX, EX, ZDUM
218 * ..
219 * .. External Functions ..
220  LOGICAL LSAME
221  INTEGER IDAMAX
222  DOUBLE PRECISION DLAMCH
223  EXTERNAL lsame, idamax, dlamch
224 * ..
225 * .. External Subroutines ..
226  EXTERNAL xerbla, zaxpy, zpttrs
227 * ..
228 * .. Intrinsic Functions ..
229  INTRINSIC abs, dble, dcmplx, dconjg, dimag, max
230 * ..
231 * .. Statement Functions ..
232  DOUBLE PRECISION CABS1
233 * ..
234 * .. Statement Function definitions ..
235  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
236 * ..
237 * .. Executable Statements ..
238 *
239 * Test the input parameters.
240 *
241  info = 0
242  upper = lsame( uplo, 'U' )
243  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
244  info = -1
245  ELSE IF( n.LT.0 ) THEN
246  info = -2
247  ELSE IF( nrhs.LT.0 ) THEN
248  info = -3
249  ELSE IF( ldb.LT.max( 1, n ) ) THEN
250  info = -9
251  ELSE IF( ldx.LT.max( 1, n ) ) THEN
252  info = -11
253  END IF
254  IF( info.NE.0 ) THEN
255  CALL xerbla( 'ZPTRFS', -info )
256  RETURN
257  END IF
258 *
259 * Quick return if possible
260 *
261  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
262  DO 10 j = 1, nrhs
263  ferr( j ) = zero
264  berr( j ) = zero
265  10 CONTINUE
266  RETURN
267  END IF
268 *
269 * NZ = maximum number of nonzero elements in each row of A, plus 1
270 *
271  nz = 4
272  eps = dlamch( 'Epsilon' )
273  safmin = dlamch( 'Safe minimum' )
274  safe1 = nz*safmin
275  safe2 = safe1 / eps
276 *
277 * Do for each right hand side
278 *
279  DO 100 j = 1, nrhs
280 *
281  count = 1
282  lstres = three
283  20 CONTINUE
284 *
285 * Loop until stopping criterion is satisfied.
286 *
287 * Compute residual R = B - A * X. Also compute
288 * abs(A)*abs(x) + abs(b) for use in the backward error bound.
289 *
290  IF( upper ) THEN
291  IF( n.EQ.1 ) THEN
292  bi = b( 1, j )
293  dx = d( 1 )*x( 1, j )
294  work( 1 ) = bi - dx
295  rwork( 1 ) = cabs1( bi ) + cabs1( dx )
296  ELSE
297  bi = b( 1, j )
298  dx = d( 1 )*x( 1, j )
299  ex = e( 1 )*x( 2, j )
300  work( 1 ) = bi - dx - ex
301  rwork( 1 ) = cabs1( bi ) + cabs1( dx ) +
302  $ cabs1( e( 1 ) )*cabs1( x( 2, j ) )
303  DO 30 i = 2, n - 1
304  bi = b( i, j )
305  cx = dconjg( e( i-1 ) )*x( i-1, j )
306  dx = d( i )*x( i, j )
307  ex = e( i )*x( i+1, j )
308  work( i ) = bi - cx - dx - ex
309  rwork( i ) = cabs1( bi ) +
310  $ cabs1( e( i-1 ) )*cabs1( x( i-1, j ) ) +
311  $ cabs1( dx ) + cabs1( e( i ) )*
312  $ cabs1( x( i+1, j ) )
313  30 CONTINUE
314  bi = b( n, j )
315  cx = dconjg( e( n-1 ) )*x( n-1, j )
316  dx = d( n )*x( n, j )
317  work( n ) = bi - cx - dx
318  rwork( n ) = cabs1( bi ) + cabs1( e( n-1 ) )*
319  $ cabs1( x( n-1, j ) ) + cabs1( dx )
320  END IF
321  ELSE
322  IF( n.EQ.1 ) THEN
323  bi = b( 1, j )
324  dx = d( 1 )*x( 1, j )
325  work( 1 ) = bi - dx
326  rwork( 1 ) = cabs1( bi ) + cabs1( dx )
327  ELSE
328  bi = b( 1, j )
329  dx = d( 1 )*x( 1, j )
330  ex = dconjg( e( 1 ) )*x( 2, j )
331  work( 1 ) = bi - dx - ex
332  rwork( 1 ) = cabs1( bi ) + cabs1( dx ) +
333  $ cabs1( e( 1 ) )*cabs1( x( 2, j ) )
334  DO 40 i = 2, n - 1
335  bi = b( i, j )
336  cx = e( i-1 )*x( i-1, j )
337  dx = d( i )*x( i, j )
338  ex = dconjg( e( i ) )*x( i+1, j )
339  work( i ) = bi - cx - dx - ex
340  rwork( i ) = cabs1( bi ) +
341  $ cabs1( e( i-1 ) )*cabs1( x( i-1, j ) ) +
342  $ cabs1( dx ) + cabs1( e( i ) )*
343  $ cabs1( x( i+1, j ) )
344  40 CONTINUE
345  bi = b( n, j )
346  cx = e( n-1 )*x( n-1, j )
347  dx = d( n )*x( n, j )
348  work( n ) = bi - cx - dx
349  rwork( n ) = cabs1( bi ) + cabs1( e( n-1 ) )*
350  $ cabs1( x( n-1, j ) ) + cabs1( dx )
351  END IF
352  END IF
353 *
354 * Compute componentwise relative backward error from formula
355 *
356 * max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
357 *
358 * where abs(Z) is the componentwise absolute value of the matrix
359 * or vector Z. If the i-th component of the denominator is less
360 * than SAFE2, then SAFE1 is added to the i-th components of the
361 * numerator and denominator before dividing.
362 *
363  s = zero
364  DO 50 i = 1, n
365  IF( rwork( i ).GT.safe2 ) THEN
366  s = max( s, cabs1( work( i ) ) / rwork( i ) )
367  ELSE
368  s = max( s, ( cabs1( work( i ) )+safe1 ) /
369  $ ( rwork( i )+safe1 ) )
370  END IF
371  50 CONTINUE
372  berr( j ) = s
373 *
374 * Test stopping criterion. Continue iterating if
375 * 1) The residual BERR(J) is larger than machine epsilon, and
376 * 2) BERR(J) decreased by at least a factor of 2 during the
377 * last iteration, and
378 * 3) At most ITMAX iterations tried.
379 *
380  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
381  $ count.LE.itmax ) THEN
382 *
383 * Update solution and try again.
384 *
385  CALL zpttrs( uplo, n, 1, df, ef, work, n, info )
386  CALL zaxpy( n, dcmplx( one ), work, 1, x( 1, j ), 1 )
387  lstres = berr( j )
388  count = count + 1
389  GO TO 20
390  END IF
391 *
392 * Bound error from formula
393 *
394 * norm(X - XTRUE) / norm(X) .le. FERR =
395 * norm( abs(inv(A))*
396 * ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
397 *
398 * where
399 * norm(Z) is the magnitude of the largest component of Z
400 * inv(A) is the inverse of A
401 * abs(Z) is the componentwise absolute value of the matrix or
402 * vector Z
403 * NZ is the maximum number of nonzeros in any row of A, plus 1
404 * EPS is machine epsilon
405 *
406 * The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
407 * is incremented by SAFE1 if the i-th component of
408 * abs(A)*abs(X) + abs(B) is less than SAFE2.
409 *
410  DO 60 i = 1, n
411  IF( rwork( i ).GT.safe2 ) THEN
412  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
413  ELSE
414  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
415  $ safe1
416  END IF
417  60 CONTINUE
418  ix = idamax( n, rwork, 1 )
419  ferr( j ) = rwork( ix )
420 *
421 * Estimate the norm of inv(A).
422 *
423 * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
424 *
425 * m(i,j) = abs(A(i,j)), i = j,
426 * m(i,j) = -abs(A(i,j)), i .ne. j,
427 *
428 * and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**H.
429 *
430 * Solve M(L) * x = e.
431 *
432  rwork( 1 ) = one
433  DO 70 i = 2, n
434  rwork( i ) = one + rwork( i-1 )*abs( ef( i-1 ) )
435  70 CONTINUE
436 *
437 * Solve D * M(L)**H * x = b.
438 *
439  rwork( n ) = rwork( n ) / df( n )
440  DO 80 i = n - 1, 1, -1
441  rwork( i ) = rwork( i ) / df( i ) +
442  $ rwork( i+1 )*abs( ef( i ) )
443  80 CONTINUE
444 *
445 * Compute norm(inv(A)) = max(x(i)), 1<=i<=n.
446 *
447  ix = idamax( n, rwork, 1 )
448  ferr( j ) = ferr( j )*abs( rwork( ix ) )
449 *
450 * Normalize error.
451 *
452  lstres = zero
453  DO 90 i = 1, n
454  lstres = max( lstres, abs( x( i, j ) ) )
455  90 CONTINUE
456  IF( lstres.NE.zero )
457  $ ferr( j ) = ferr( j ) / lstres
458 *
459  100 CONTINUE
460 *
461  RETURN
462 *
463 * End of ZPTRFS
464 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:88
subroutine zpttrs(UPLO, N, NRHS, D, E, B, LDB, INFO)
ZPTTRS
Definition: zpttrs.f:121
Here is the call graph for this function:
Here is the caller graph for this function: