LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ csprfs()

subroutine csprfs ( character  UPLO,
integer  N,
integer  NRHS,
complex, dimension( * )  AP,
complex, dimension( * )  AFP,
integer, dimension( * )  IPIV,
complex, dimension( ldb, * )  B,
integer  LDB,
complex, dimension( ldx, * )  X,
integer  LDX,
real, dimension( * )  FERR,
real, dimension( * )  BERR,
complex, dimension( * )  WORK,
real, dimension( * )  RWORK,
integer  INFO 
)

CSPRFS

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

Purpose:
 CSPRFS 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 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 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 CSPTRF, 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 CSPTRF.
[in]B
          B is COMPLEX 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 array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by CSPTRS.
          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 COMPLEX array, dimension (2*N)
[out]RWORK
          RWORK is REAL 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 178 of file csprfs.f.

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