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

◆ 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dpttrs
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: