LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
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.`

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