LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2011

Definition at line 173 of file dpprfs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: