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

◆ dpprfs()

subroutine dpprfs ( character  uplo,
integer  n,
integer  nrhs,
double precision, dimension( * )  ap,
double precision, dimension( * )  afp,
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 
)

DPPRFS

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

Purpose:
 DPPRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is symmetric positive definite
 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)*(2n-j)/2) = A(i,j) for j<=i<=n.
[in]AFP
          AFP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
          The triangular factor U or L from the Cholesky factorization
          A = U**T*U or A = L*L**T, as computed by DPPTRF/ZPPTRF,
          packed columnwise in a linear array in the same format as A
          (see AP).
[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 DPPTRS.
          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 169 of file dpprfs.f.

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