LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dptrfs()

subroutine dptrfs ( integer  N,
integer  NRHS,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( * )  DF,
double precision, dimension( * )  EF,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( ldx, * )  X,
integer  LDX,
double precision, dimension( * )  FERR,
double precision, dimension( * )  BERR,
double precision, dimension( * )  WORK,
integer  INFO 
)

DPTRFS

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

Purpose:
 DPTRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is symmetric positive definite
 and tridiagonal, and provides error bounds and backward error
 estimates for the solution.
Parameters
[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 diagonal elements of the tridiagonal matrix A.
[in]E
          E is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) subdiagonal elements of the tridiagonal matrix A.
[in]DF
          DF is DOUBLE PRECISION array, dimension (N)
          The n diagonal elements of the diagonal matrix D from the
          factorization computed by DPTTRF.
[in]EF
          EF is DOUBLE PRECISION array, dimension (N-1)
          The (n-1) subdiagonal elements of the unit bidiagonal factor
          L from the factorization computed by DPTTRF.
[in]B
          B is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by DPTTRS.
          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 DOUBLE PRECISION array, dimension (2*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 161 of file dptrfs.f.

163 *
164 * -- LAPACK computational routine --
165 * -- LAPACK is a software package provided by Univ. of Tennessee, --
166 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
167 *
168 * .. Scalar Arguments ..
169  INTEGER INFO, LDB, LDX, N, NRHS
170 * ..
171 * .. Array Arguments ..
172  DOUBLE PRECISION B( LDB, * ), BERR( * ), D( * ), DF( * ),
173  $ E( * ), EF( * ), FERR( * ), WORK( * ),
174  $ X( LDX, * )
175 * ..
176 *
177 * =====================================================================
178 *
179 * .. Parameters ..
180  INTEGER ITMAX
181  parameter( itmax = 5 )
182  DOUBLE PRECISION ZERO
183  parameter( zero = 0.0d+0 )
184  DOUBLE PRECISION ONE
185  parameter( one = 1.0d+0 )
186  DOUBLE PRECISION TWO
187  parameter( two = 2.0d+0 )
188  DOUBLE PRECISION THREE
189  parameter( three = 3.0d+0 )
190 * ..
191 * .. Local Scalars ..
192  INTEGER COUNT, I, IX, J, NZ
193  DOUBLE PRECISION BI, CX, DX, EPS, EX, LSTRES, S, SAFE1, SAFE2,
194  $ SAFMIN
195 * ..
196 * .. External Subroutines ..
197  EXTERNAL daxpy, dpttrs, xerbla
198 * ..
199 * .. Intrinsic Functions ..
200  INTRINSIC abs, max
201 * ..
202 * .. External Functions ..
203  INTEGER IDAMAX
204  DOUBLE PRECISION DLAMCH
205  EXTERNAL idamax, dlamch
206 * ..
207 * .. Executable Statements ..
208 *
209 * Test the input parameters.
210 *
211  info = 0
212  IF( n.LT.0 ) THEN
213  info = -1
214  ELSE IF( nrhs.LT.0 ) THEN
215  info = -2
216  ELSE IF( ldb.LT.max( 1, n ) ) THEN
217  info = -8
218  ELSE IF( ldx.LT.max( 1, n ) ) THEN
219  info = -10
220  END IF
221  IF( info.NE.0 ) THEN
222  CALL xerbla( 'DPTRFS', -info )
223  RETURN
224  END IF
225 *
226 * Quick return if possible
227 *
228  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
229  DO 10 j = 1, nrhs
230  ferr( j ) = zero
231  berr( j ) = zero
232  10 CONTINUE
233  RETURN
234  END IF
235 *
236 * NZ = maximum number of nonzero elements in each row of A, plus 1
237 *
238  nz = 4
239  eps = dlamch( 'Epsilon' )
240  safmin = dlamch( 'Safe minimum' )
241  safe1 = nz*safmin
242  safe2 = safe1 / eps
243 *
244 * Do for each right hand side
245 *
246  DO 90 j = 1, nrhs
247 *
248  count = 1
249  lstres = three
250  20 CONTINUE
251 *
252 * Loop until stopping criterion is satisfied.
253 *
254 * Compute residual R = B - A * X. Also compute
255 * abs(A)*abs(x) + abs(b) for use in the backward error bound.
256 *
257  IF( n.EQ.1 ) THEN
258  bi = b( 1, j )
259  dx = d( 1 )*x( 1, j )
260  work( n+1 ) = bi - dx
261  work( 1 ) = abs( bi ) + abs( dx )
262  ELSE
263  bi = b( 1, j )
264  dx = d( 1 )*x( 1, j )
265  ex = e( 1 )*x( 2, j )
266  work( n+1 ) = bi - dx - ex
267  work( 1 ) = abs( bi ) + abs( dx ) + abs( ex )
268  DO 30 i = 2, n - 1
269  bi = b( i, j )
270  cx = e( i-1 )*x( i-1, j )
271  dx = d( i )*x( i, j )
272  ex = e( i )*x( i+1, j )
273  work( n+i ) = bi - cx - dx - ex
274  work( i ) = abs( bi ) + abs( cx ) + abs( dx ) + abs( ex )
275  30 CONTINUE
276  bi = b( n, j )
277  cx = e( n-1 )*x( n-1, j )
278  dx = d( n )*x( n, j )
279  work( n+n ) = bi - cx - dx
280  work( n ) = abs( bi ) + abs( cx ) + abs( dx )
281  END IF
282 *
283 * Compute componentwise relative backward error from formula
284 *
285 * max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
286 *
287 * where abs(Z) is the componentwise absolute value of the matrix
288 * or vector Z. If the i-th component of the denominator is less
289 * than SAFE2, then SAFE1 is added to the i-th components of the
290 * numerator and denominator before dividing.
291 *
292  s = zero
293  DO 40 i = 1, n
294  IF( work( i ).GT.safe2 ) THEN
295  s = max( s, abs( work( n+i ) ) / work( i ) )
296  ELSE
297  s = max( s, ( abs( work( n+i ) )+safe1 ) /
298  $ ( work( i )+safe1 ) )
299  END IF
300  40 CONTINUE
301  berr( j ) = s
302 *
303 * Test stopping criterion. Continue iterating if
304 * 1) The residual BERR(J) is larger than machine epsilon, and
305 * 2) BERR(J) decreased by at least a factor of 2 during the
306 * last iteration, and
307 * 3) At most ITMAX iterations tried.
308 *
309  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
310  $ count.LE.itmax ) THEN
311 *
312 * Update solution and try again.
313 *
314  CALL dpttrs( n, 1, df, ef, work( n+1 ), n, info )
315  CALL daxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
316  lstres = berr( j )
317  count = count + 1
318  GO TO 20
319  END IF
320 *
321 * Bound error from formula
322 *
323 * norm(X - XTRUE) / norm(X) .le. FERR =
324 * norm( abs(inv(A))*
325 * ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
326 *
327 * where
328 * norm(Z) is the magnitude of the largest component of Z
329 * inv(A) is the inverse of A
330 * abs(Z) is the componentwise absolute value of the matrix or
331 * vector Z
332 * NZ is the maximum number of nonzeros in any row of A, plus 1
333 * EPS is machine epsilon
334 *
335 * The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
336 * is incremented by SAFE1 if the i-th component of
337 * abs(A)*abs(X) + abs(B) is less than SAFE2.
338 *
339  DO 50 i = 1, n
340  IF( work( i ).GT.safe2 ) THEN
341  work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
342  ELSE
343  work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
344  END IF
345  50 CONTINUE
346  ix = idamax( n, work, 1 )
347  ferr( j ) = work( ix )
348 *
349 * Estimate the norm of inv(A).
350 *
351 * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
352 *
353 * m(i,j) = abs(A(i,j)), i = j,
354 * m(i,j) = -abs(A(i,j)), i .ne. j,
355 *
356 * and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**T.
357 *
358 * Solve M(L) * x = e.
359 *
360  work( 1 ) = one
361  DO 60 i = 2, n
362  work( i ) = one + work( i-1 )*abs( ef( i-1 ) )
363  60 CONTINUE
364 *
365 * Solve D * M(L)**T * x = b.
366 *
367  work( n ) = work( n ) / df( n )
368  DO 70 i = n - 1, 1, -1
369  work( i ) = work( i ) / df( i ) + work( i+1 )*abs( ef( i ) )
370  70 CONTINUE
371 *
372 * Compute norm(inv(A)) = max(x(i)), 1<=i<=n.
373 *
374  ix = idamax( n, work, 1 )
375  ferr( j ) = ferr( j )*abs( work( ix ) )
376 *
377 * Normalize error.
378 *
379  lstres = zero
380  DO 80 i = 1, n
381  lstres = max( lstres, abs( x( i, j ) ) )
382  80 CONTINUE
383  IF( lstres.NE.zero )
384  $ ferr( j ) = ferr( j ) / lstres
385 *
386  90 CONTINUE
387 *
388  RETURN
389 *
390 * End of DPTRFS
391 *
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
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dpttrs(N, NRHS, D, E, B, LDB, INFO)
DPTTRS
Definition: dpttrs.f:109
Here is the call graph for this function:
Here is the caller graph for this function: