LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ spprfs()

subroutine spprfs ( character  UPLO,
integer  N,
integer  NRHS,
real, dimension( * )  AP,
real, dimension( * )  AFP,
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 
)

SPPRFS

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

Purpose:
 SPPRFS 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 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)*(2n-j)/2) = A(i,j) for j<=i<=n.
[in]AFP
          AFP is REAL 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 SPPTRF/CPPTRF,
          packed columnwise in a linear array in the same format as A
          (see AP).
[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 SPPTRS.
          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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 169 of file spprfs.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  REAL AFP( * ), AP( * ), B( LDB, * ), BERR( * ),
183  $ FERR( * ), WORK( * ), X( LDX, * )
184 * ..
185 *
186 * =====================================================================
187 *
188 * .. Parameters ..
189  INTEGER ITMAX
190  parameter( itmax = 5 )
191  REAL ZERO
192  parameter( zero = 0.0e+0 )
193  REAL ONE
194  parameter( one = 1.0e+0 )
195  REAL TWO
196  parameter( two = 2.0e+0 )
197  REAL THREE
198  parameter( three = 3.0e+0 )
199 * ..
200 * .. Local Scalars ..
201  LOGICAL UPPER
202  INTEGER COUNT, I, IK, J, K, KASE, KK, NZ
203  REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
204 * ..
205 * .. Local Arrays ..
206  INTEGER ISAVE( 3 )
207 * ..
208 * .. External Subroutines ..
209  EXTERNAL saxpy, scopy, slacn2, spptrs, sspmv, xerbla
210 * ..
211 * .. Intrinsic Functions ..
212  INTRINSIC abs, max
213 * ..
214 * .. External Functions ..
215  LOGICAL LSAME
216  REAL SLAMCH
217  EXTERNAL lsame, slamch
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( 'SPPRFS', -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 = slamch( 'Epsilon' )
255  safmin = slamch( '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 scopy( n, b( 1, j ), 1, work( n+1 ), 1 )
272  CALL sspmv( 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 spptrs( uplo, n, 1, afp, work( n+1 ), n, info )
342  CALL saxpy( 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 SLACN2 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 slacn2( 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 spptrs( 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 spptrs( 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 SPPRFS
417 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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:136
subroutine spptrs(UPLO, N, NRHS, AP, B, LDB, INFO)
SPPTRS
Definition: spptrs.f:108
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89
subroutine sspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
SSPMV
Definition: sspmv.f:147
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: