LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zsprfs ( character  UPLO,
integer  N,
integer  NRHS,
complex*16, dimension( * )  AP,
complex*16, dimension( * )  AFP,
integer, dimension( * )  IPIV,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldx, * )  X,
integer  LDX,
double precision, dimension( * )  FERR,
double precision, dimension( * )  BERR,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZSPRFS

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

Purpose: 
 ZSPRFS 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 COMPLEX*16 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 COMPLEX*16 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 ZSPTRF, 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 ZSPTRF.
[in]B
          B is COMPLEX*16 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 COMPLEX*16 array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by ZSPTRS.
          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 COMPLEX*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION 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 182 of file zsprfs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: