LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 subroutine sptrfs ( integer N, integer NRHS, real, dimension( * ) D, real, dimension( * ) E, real, dimension( * ) DF, real, dimension( * ) EF, real, dimension( ldb, * ) B, integer LDB, real, dimension( ldx, * ) X, integer LDX, real, dimension( * ) FERR, real, dimension( * ) BERR, real, dimension( * ) WORK, integer INFO )

SPTRFS

Purpose:
``` SPTRFS 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 REAL array, dimension (N) The n diagonal elements of the tridiagonal matrix A.``` [in] E ``` E is REAL array, dimension (N-1) The (n-1) subdiagonal elements of the tridiagonal matrix A.``` [in] DF ``` DF is REAL array, dimension (N) The n diagonal elements of the diagonal matrix D from the factorization computed by SPTTRF.``` [in] EF ``` EF is REAL array, dimension (N-1) The (n-1) subdiagonal elements of the unit bidiagonal factor L from the factorization computed by SPTTRF.``` [in] B ``` B is REAL 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 REAL array, dimension (LDX,NRHS) On entry, the solution matrix X, as computed by SPTTRS. 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 REAL 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 REAL 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 REAL 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.`
Date
September 2012

Definition at line 165 of file sptrfs.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  REAL 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  REAL zero
186  parameter ( zero = 0.0e+0 )
187  REAL one
188  parameter ( one = 1.0e+0 )
189  REAL two
190  parameter ( two = 2.0e+0 )
191  REAL three
192  parameter ( three = 3.0e+0 )
193 * ..
194 * .. Local Scalars ..
195  INTEGER count, i, ix, j, nz
196  REAL bi, cx, dx, eps, ex, lstres, s, safe1, safe2,
197  \$ safmin
198 * ..
199 * .. External Subroutines ..
200  EXTERNAL saxpy, spttrs, xerbla
201 * ..
202 * .. Intrinsic Functions ..
203  INTRINSIC abs, max
204 * ..
205 * .. External Functions ..
206  INTEGER isamax
207  REAL slamch
208  EXTERNAL isamax, slamch
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( 'SPTRFS', -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 = slamch( 'Epsilon' )
243  safmin = slamch( '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 spttrs( n, 1, df, ef, work( n+1 ), n, info )
318  CALL saxpy( 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 = isamax( 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 = isamax( 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 SPTRFS
394 *
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:54
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine spttrs(N, NRHS, D, E, B, LDB, INFO)
SPTTRS
Definition: spttrs.f:111

Here is the call graph for this function:

Here is the caller graph for this function:
{"ft_lang_label":"__label__en","ft_lang_prob":0.71967936,"math_prob":0.9975309,"size":2544,"snap":"2019-26-2019-30","text_gpt3_token_len":737,"char_repetition_ratio":0.17125984,"word_repetition_ratio":0.13333334,"special_character_ratio":0.2841981,"punctuation_ratio":0.15204678,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9990332,"pos_list":[0,1,2],"im_url_duplicate_count":[null,null,null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2019-07-21T07:20:07Z\",\"WARC-Record-ID\":\"<urn:uuid:cbf7fde2-28e8-4509-8d65-5170c1731eff>\",\"Content-Length\":\"50279\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:3235d00f-2e57-446c-b854-c5be3baa48ee>\",\"WARC-Concurrent-To\":\"<urn:uuid:ea84f626-7af6-4cfe-9c45-4b28d9cccc7f>\",\"WARC-IP-Address\":\"160.36.131.221\",\"WARC-Target-URI\":\"http://www.netlib.org/lapack/explore-html-3.6.1/d1/ddc/group__real_p_tcomputational_ga50fb1dfb45818407f7748adc9f804f58.html\",\"WARC-Payload-Digest\":\"sha1:RHBPRT2KKQPJ4H7LQQPGBPPX6W2E2MKX\",\"WARC-Block-Digest\":\"sha1:EQ7IDCPHAGZGL3JEWYR2CCONYL2UYP4X\",\"WARC-Identified-Payload-Type\":\"application/xhtml+xml\",\"warc_filename\":\"/cc_download/warc_2019/CC-MAIN-2019-30/CC-MAIN-2019-30_segments_1563195526931.25_warc_CC-MAIN-20190721061720-20190721083720-00537.warc.gz\"}"}