LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
September 2012

Definition at line 165 of file dptrfs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: