 LAPACK  3.8.0 LAPACK: Linear Algebra PACKage

## ◆ ssprfs()

 subroutine ssprfs ( character UPLO, integer N, integer NRHS, real, dimension( * ) AP, real, dimension( * ) AFP, integer, dimension( * ) IPIV, real, dimension( ldb, * ) B, integer LDB, real, dimension( ldx, * ) X, integer LDX, real, dimension( * ) FERR, real, dimension( * ) BERR, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO )

SSPRFS

Purpose:
``` SSPRFS improves the computed solution to a system of linear
equations when the coefficient matrix is symmetric indefinite
and packed, and provides error bounds and backward error estimates
for the solution.```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored.``` [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 matrices B and X. NRHS >= 0.``` [in] AP ``` AP is REAL array, dimension (N*(N+1)/2) The upper or lower triangle of the symmetric matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows: if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.``` [in] AFP ``` AFP is REAL array, dimension (N*(N+1)/2) The factored form of the matrix A. AFP contains the block diagonal matrix D and the multipliers used to obtain the factor U or L from the factorization A = U*D*U**T or A = L*D*L**T as computed by SSPTRF, stored as a packed triangular matrix.``` [in] IPIV ``` IPIV is INTEGER array, dimension (N) Details of the interchanges and the block structure of D as determined by SSPTRF.``` [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 SSPTRS. 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 estimated 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). The estimate is as reliable as the estimate for RCOND, and is almost always a slight overestimate of the true error.``` [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 (3*N)` [out] IWORK ` IWORK is INTEGER array, dimension (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
December 2016

Definition at line 181 of file ssprfs.f.

181 *
182 * -- LAPACK computational routine (version 3.7.0) --
183 * -- LAPACK is a software package provided by Univ. of Tennessee, --
184 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
185 * December 2016
186 *
187 * .. Scalar Arguments ..
188  CHARACTER uplo
189  INTEGER info, ldb, ldx, n, nrhs
190 * ..
191 * .. Array Arguments ..
192  INTEGER ipiv( * ), iwork( * )
193  REAL afp( * ), ap( * ), b( ldb, * ), berr( * ),
194  \$ ferr( * ), work( * ), x( ldx, * )
195 * ..
196 *
197 * =====================================================================
198 *
199 * .. Parameters ..
200  INTEGER itmax
201  parameter( itmax = 5 )
202  REAL zero
203  parameter( zero = 0.0e+0 )
204  REAL one
205  parameter( one = 1.0e+0 )
206  REAL two
207  parameter( two = 2.0e+0 )
208  REAL three
209  parameter( three = 3.0e+0 )
210 * ..
211 * .. Local Scalars ..
212  LOGICAL upper
213  INTEGER count, i, ik, j, k, kase, kk, nz
214  REAL eps, lstres, s, safe1, safe2, safmin, xk
215 * ..
216 * .. Local Arrays ..
217  INTEGER isave( 3 )
218 * ..
219 * .. External Subroutines ..
220  EXTERNAL saxpy, scopy, slacn2, sspmv, ssptrs, xerbla
221 * ..
222 * .. Intrinsic Functions ..
223  INTRINSIC abs, max
224 * ..
225 * .. External Functions ..
226  LOGICAL lsame
227  REAL slamch
228  EXTERNAL lsame, slamch
229 * ..
230 * .. Executable Statements ..
231 *
232 * Test the input parameters.
233 *
234  info = 0
235  upper = lsame( uplo, 'U' )
236  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
237  info = -1
238  ELSE IF( n.LT.0 ) THEN
239  info = -2
240  ELSE IF( nrhs.LT.0 ) THEN
241  info = -3
242  ELSE IF( ldb.LT.max( 1, n ) ) THEN
243  info = -8
244  ELSE IF( ldx.LT.max( 1, n ) ) THEN
245  info = -10
246  END IF
247  IF( info.NE.0 ) THEN
248  CALL xerbla( 'SSPRFS', -info )
249  RETURN
250  END IF
251 *
252 * Quick return if possible
253 *
254  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
255  DO 10 j = 1, nrhs
256  ferr( j ) = zero
257  berr( j ) = zero
258  10 CONTINUE
259  RETURN
260  END IF
261 *
262 * NZ = maximum number of nonzero elements in each row of A, plus 1
263 *
264  nz = n + 1
265  eps = slamch( 'Epsilon' )
266  safmin = slamch( 'Safe minimum' )
267  safe1 = nz*safmin
268  safe2 = safe1 / eps
269 *
270 * Do for each right hand side
271 *
272  DO 140 j = 1, nrhs
273 *
274  count = 1
275  lstres = three
276  20 CONTINUE
277 *
278 * Loop until stopping criterion is satisfied.
279 *
280 * Compute residual R = B - A * X
281 *
282  CALL scopy( n, b( 1, j ), 1, work( n+1 ), 1 )
283  CALL sspmv( uplo, n, -one, ap, x( 1, j ), 1, one, work( n+1 ),
284  \$ 1 )
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  DO 30 i = 1, n
296  work( i ) = abs( b( i, j ) )
297  30 CONTINUE
298 *
299 * Compute abs(A)*abs(X) + abs(B).
300 *
301  kk = 1
302  IF( upper ) THEN
303  DO 50 k = 1, n
304  s = zero
305  xk = abs( x( k, j ) )
306  ik = kk
307  DO 40 i = 1, k - 1
308  work( i ) = work( i ) + abs( ap( ik ) )*xk
309  s = s + abs( ap( ik ) )*abs( x( i, j ) )
310  ik = ik + 1
311  40 CONTINUE
312  work( k ) = work( k ) + abs( ap( kk+k-1 ) )*xk + s
313  kk = kk + k
314  50 CONTINUE
315  ELSE
316  DO 70 k = 1, n
317  s = zero
318  xk = abs( x( k, j ) )
319  work( k ) = work( k ) + abs( ap( kk ) )*xk
320  ik = kk + 1
321  DO 60 i = k + 1, n
322  work( i ) = work( i ) + abs( ap( ik ) )*xk
323  s = s + abs( ap( ik ) )*abs( x( i, j ) )
324  ik = ik + 1
325  60 CONTINUE
326  work( k ) = work( k ) + s
327  kk = kk + ( n-k+1 )
328  70 CONTINUE
329  END IF
330  s = zero
331  DO 80 i = 1, n
332  IF( work( i ).GT.safe2 ) THEN
333  s = max( s, abs( work( n+i ) ) / work( i ) )
334  ELSE
335  s = max( s, ( abs( work( n+i ) )+safe1 ) /
336  \$ ( work( i )+safe1 ) )
337  END IF
338  80 CONTINUE
339  berr( j ) = s
340 *
341 * Test stopping criterion. Continue iterating if
342 * 1) The residual BERR(J) is larger than machine epsilon, and
343 * 2) BERR(J) decreased by at least a factor of 2 during the
344 * last iteration, and
345 * 3) At most ITMAX iterations tried.
346 *
347  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
348  \$ count.LE.itmax ) THEN
349 *
350 * Update solution and try again.
351 *
352  CALL ssptrs( uplo, n, 1, afp, ipiv, work( n+1 ), n, info )
353  CALL saxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
354  lstres = berr( j )
355  count = count + 1
356  GO TO 20
357  END IF
358 *
359 * Bound error from formula
360 *
361 * norm(X - XTRUE) / norm(X) .le. FERR =
362 * norm( abs(inv(A))*
363 * ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
364 *
365 * where
366 * norm(Z) is the magnitude of the largest component of Z
367 * inv(A) is the inverse of A
368 * abs(Z) is the componentwise absolute value of the matrix or
369 * vector Z
370 * NZ is the maximum number of nonzeros in any row of A, plus 1
371 * EPS is machine epsilon
372 *
373 * The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
374 * is incremented by SAFE1 if the i-th component of
375 * abs(A)*abs(X) + abs(B) is less than SAFE2.
376 *
377 * Use SLACN2 to estimate the infinity-norm of the matrix
378 * inv(A) * diag(W),
379 * where W = abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) )))
380 *
381  DO 90 i = 1, n
382  IF( work( i ).GT.safe2 ) THEN
383  work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
384  ELSE
385  work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
386  END IF
387  90 CONTINUE
388 *
389  kase = 0
390  100 CONTINUE
391  CALL slacn2( n, work( 2*n+1 ), work( n+1 ), iwork, ferr( j ),
392  \$ kase, isave )
393  IF( kase.NE.0 ) THEN
394  IF( kase.EQ.1 ) THEN
395 *
396 * Multiply by diag(W)*inv(A**T).
397 *
398  CALL ssptrs( uplo, n, 1, afp, ipiv, work( n+1 ), n,
399  \$ info )
400  DO 110 i = 1, n
401  work( n+i ) = work( i )*work( n+i )
402  110 CONTINUE
403  ELSE IF( kase.EQ.2 ) THEN
404 *
405 * Multiply by inv(A)*diag(W).
406 *
407  DO 120 i = 1, n
408  work( n+i ) = work( i )*work( n+i )
409  120 CONTINUE
410  CALL ssptrs( uplo, n, 1, afp, ipiv, work( n+1 ), n,
411  \$ info )
412  END IF
413  GO TO 100
414  END IF
415 *
416 * Normalize error.
417 *
418  lstres = zero
419  DO 130 i = 1, n
420  lstres = max( lstres, abs( x( i, j ) ) )
421  130 CONTINUE
422  IF( lstres.NE.zero )
423  \$ ferr( j ) = ferr( j ) / lstres
424 *
425  140 CONTINUE
426 *
427  RETURN
428 *
429 * End of SSPRFS
430 *
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:91
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine sspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
SSPMV
Definition: sspmv.f:149
subroutine slacn2(N, V, X, ISGN, EST, KASE, ISAVE)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: slacn2.f:138
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine ssptrs(UPLO, N, NRHS, AP, IPIV, B, LDB, INFO)
SSPTRS
Definition: ssptrs.f:117
Here is the call graph for this function:
Here is the caller graph for this function: