LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ dsprfs()

subroutine dsprfs ( character  UPLO,
integer  N,
integer  NRHS,
double precision, dimension( * )  AP,
double precision, dimension( * )  AFP,
integer, dimension( * )  IPIV,
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, dimension( * )  IWORK,
integer  INFO 
)

DSPRFS

Download DSPRFS + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DSPRFS 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DSPTRF, 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 DSPTRF.
[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 DSPTRS.
          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 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 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 (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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 177 of file dsprfs.f.

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