LAPACK 3.12.1
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 159 of file dptrfs.f.

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