LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dspmv
subroutine dspmv(uplo, n, alpha, ap, x, incx, beta, y, incy)
DSPMV
Definition dspmv.f:147
dsptrs
subroutine dsptrs(uplo, n, nrhs, ap, ipiv, b, ldb, info)
DSPTRS
Definition dsptrs.f:115
dlacn2
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
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the call graph for this function:
Here is the caller graph for this function: