LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ctprfs ( character  UPLO,
character  TRANS,
character  DIAG,
integer  N,
integer  NRHS,
complex, dimension( * )  AP,
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 
)

CTPRFS

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

Purpose: 
 CTPRFS provides error bounds and backward error estimates for the
 solution to a system of linear equations with a triangular packed
 coefficient matrix.

 The solution matrix X must be computed by CTPTRS or some other
 means before entering this routine.  CTPRFS does not do iterative
 refinement because doing so cannot improve the backward error.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  A is upper triangular;
          = 'L':  A is lower triangular.
[in]TRANS
          TRANS is CHARACTER*1
          Specifies the form of the system of equations:
          = 'N':  A * X = B     (No transpose)
          = 'T':  A**T * X = B  (Transpose)
          = 'C':  A**H * X = B  (Conjugate transpose)
[in]DIAG
          DIAG is CHARACTER*1
          = 'N':  A is non-unit triangular;
          = 'U':  A is unit triangular.
[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 triangular 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.
          If DIAG = 'U', the diagonal elements of A are not referenced
          and are assumed to be 1.
[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]X
          X is COMPLEX array, dimension (LDX,NRHS)
          The 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
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 176 of file ctprfs.f.

176 *
177 * -- LAPACK computational routine (version 3.4.0) --
178 * -- LAPACK is a software package provided by Univ. of Tennessee, --
179 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180 * November 2011
181 *
182 * .. Scalar Arguments ..
183  CHARACTER diag, trans, uplo
184  INTEGER info, ldb, ldx, n, nrhs
185 * ..
186 * .. Array Arguments ..
187  REAL berr( * ), ferr( * ), rwork( * )
188  COMPLEX ap( * ), b( ldb, * ), work( * ), x( ldx, * )
189 * ..
190 *
191 * =====================================================================
192 *
193 * .. Parameters ..
194  REAL zero
195  parameter ( zero = 0.0e+0 )
196  COMPLEX one
197  parameter ( one = ( 1.0e+0, 0.0e+0 ) )
198 * ..
199 * .. Local Scalars ..
200  LOGICAL notran, nounit, upper
201  CHARACTER transn, transt
202  INTEGER i, j, k, kase, kc, nz
203  REAL eps, lstres, s, safe1, safe2, safmin, xk
204  COMPLEX zdum
205 * ..
206 * .. Local Arrays ..
207  INTEGER isave( 3 )
208 * ..
209 * .. External Subroutines ..
210  EXTERNAL caxpy, ccopy, clacn2, ctpmv, ctpsv, xerbla
211 * ..
212 * .. Intrinsic Functions ..
213  INTRINSIC abs, aimag, max, real
214 * ..
215 * .. External Functions ..
216  LOGICAL lsame
217  REAL slamch
218  EXTERNAL lsame, slamch
219 * ..
220 * .. Statement Functions ..
221  REAL cabs1
222 * ..
223 * .. Statement Function definitions ..
224  cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( aimag( zdum ) )
225 * ..
226 * .. Executable Statements ..
227 *
228 * Test the input parameters.
229 *
230  info = 0
231  upper = lsame( uplo, 'U' )
232  notran = lsame( trans, 'N' )
233  nounit = lsame( diag, 'N' )
234 *
235  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
236  info = -1
237  ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
238  $ lsame( trans, 'C' ) ) THEN
239  info = -2
240  ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
241  info = -3
242  ELSE IF( n.LT.0 ) THEN
243  info = -4
244  ELSE IF( nrhs.LT.0 ) THEN
245  info = -5
246  ELSE IF( ldb.LT.max( 1, n ) ) THEN
247  info = -8
248  ELSE IF( ldx.LT.max( 1, n ) ) THEN
249  info = -10
250  END IF
251  IF( info.NE.0 ) THEN
252  CALL xerbla( 'CTPRFS', -info )
253  RETURN
254  END IF
255 *
256 * Quick return if possible
257 *
258  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
259  DO 10 j = 1, nrhs
260  ferr( j ) = zero
261  berr( j ) = zero
262  10 CONTINUE
263  RETURN
264  END IF
265 *
266  IF( notran ) THEN
267  transn = 'N'
268  transt = 'C'
269  ELSE
270  transn = 'C'
271  transt = 'N'
272  END IF
273 *
274 * NZ = maximum number of nonzero elements in each row of A, plus 1
275 *
276  nz = n + 1
277  eps = slamch( 'Epsilon' )
278  safmin = slamch( 'Safe minimum' )
279  safe1 = nz*safmin
280  safe2 = safe1 / eps
281 *
282 * Do for each right hand side
283 *
284  DO 250 j = 1, nrhs
285 *
286 * Compute residual R = B - op(A) * X,
287 * where op(A) = A, A**T, or A**H, depending on TRANS.
288 *
289  CALL ccopy( n, x( 1, j ), 1, work, 1 )
290  CALL ctpmv( uplo, trans, diag, n, ap, work, 1 )
291  CALL caxpy( n, -one, b( 1, j ), 1, work, 1 )
292 *
293 * Compute componentwise relative backward error from formula
294 *
295 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
296 *
297 * where abs(Z) is the componentwise absolute value of the matrix
298 * or vector Z. If the i-th component of the denominator is less
299 * than SAFE2, then SAFE1 is added to the i-th components of the
300 * numerator and denominator before dividing.
301 *
302  DO 20 i = 1, n
303  rwork( i ) = cabs1( b( i, j ) )
304  20 CONTINUE
305 *
306  IF( notran ) THEN
307 *
308 * Compute abs(A)*abs(X) + abs(B).
309 *
310  IF( upper ) THEN
311  kc = 1
312  IF( nounit ) THEN
313  DO 40 k = 1, n
314  xk = cabs1( x( k, j ) )
315  DO 30 i = 1, k
316  rwork( i ) = rwork( i ) +
317  $ cabs1( ap( kc+i-1 ) )*xk
318  30 CONTINUE
319  kc = kc + k
320  40 CONTINUE
321  ELSE
322  DO 60 k = 1, n
323  xk = cabs1( x( k, j ) )
324  DO 50 i = 1, k - 1
325  rwork( i ) = rwork( i ) +
326  $ cabs1( ap( kc+i-1 ) )*xk
327  50 CONTINUE
328  rwork( k ) = rwork( k ) + xk
329  kc = kc + k
330  60 CONTINUE
331  END IF
332  ELSE
333  kc = 1
334  IF( nounit ) THEN
335  DO 80 k = 1, n
336  xk = cabs1( x( k, j ) )
337  DO 70 i = k, n
338  rwork( i ) = rwork( i ) +
339  $ cabs1( ap( kc+i-k ) )*xk
340  70 CONTINUE
341  kc = kc + n - k + 1
342  80 CONTINUE
343  ELSE
344  DO 100 k = 1, n
345  xk = cabs1( x( k, j ) )
346  DO 90 i = k + 1, n
347  rwork( i ) = rwork( i ) +
348  $ cabs1( ap( kc+i-k ) )*xk
349  90 CONTINUE
350  rwork( k ) = rwork( k ) + xk
351  kc = kc + n - k + 1
352  100 CONTINUE
353  END IF
354  END IF
355  ELSE
356 *
357 * Compute abs(A**H)*abs(X) + abs(B).
358 *
359  IF( upper ) THEN
360  kc = 1
361  IF( nounit ) THEN
362  DO 120 k = 1, n
363  s = zero
364  DO 110 i = 1, k
365  s = s + cabs1( ap( kc+i-1 ) )*cabs1( x( i, j ) )
366  110 CONTINUE
367  rwork( k ) = rwork( k ) + s
368  kc = kc + k
369  120 CONTINUE
370  ELSE
371  DO 140 k = 1, n
372  s = cabs1( x( k, j ) )
373  DO 130 i = 1, k - 1
374  s = s + cabs1( ap( kc+i-1 ) )*cabs1( x( i, j ) )
375  130 CONTINUE
376  rwork( k ) = rwork( k ) + s
377  kc = kc + k
378  140 CONTINUE
379  END IF
380  ELSE
381  kc = 1
382  IF( nounit ) THEN
383  DO 160 k = 1, n
384  s = zero
385  DO 150 i = k, n
386  s = s + cabs1( ap( kc+i-k ) )*cabs1( x( i, j ) )
387  150 CONTINUE
388  rwork( k ) = rwork( k ) + s
389  kc = kc + n - k + 1
390  160 CONTINUE
391  ELSE
392  DO 180 k = 1, n
393  s = cabs1( x( k, j ) )
394  DO 170 i = k + 1, n
395  s = s + cabs1( ap( kc+i-k ) )*cabs1( x( i, j ) )
396  170 CONTINUE
397  rwork( k ) = rwork( k ) + s
398  kc = kc + n - k + 1
399  180 CONTINUE
400  END IF
401  END IF
402  END IF
403  s = zero
404  DO 190 i = 1, n
405  IF( rwork( i ).GT.safe2 ) THEN
406  s = max( s, cabs1( work( i ) ) / rwork( i ) )
407  ELSE
408  s = max( s, ( cabs1( work( i ) )+safe1 ) /
409  $ ( rwork( i )+safe1 ) )
410  END IF
411  190 CONTINUE
412  berr( j ) = s
413 *
414 * Bound error from formula
415 *
416 * norm(X - XTRUE) / norm(X) .le. FERR =
417 * norm( abs(inv(op(A)))*
418 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
419 *
420 * where
421 * norm(Z) is the magnitude of the largest component of Z
422 * inv(op(A)) is the inverse of op(A)
423 * abs(Z) is the componentwise absolute value of the matrix or
424 * vector Z
425 * NZ is the maximum number of nonzeros in any row of A, plus 1
426 * EPS is machine epsilon
427 *
428 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
429 * is incremented by SAFE1 if the i-th component of
430 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
431 *
432 * Use CLACN2 to estimate the infinity-norm of the matrix
433 * inv(op(A)) * diag(W),
434 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
435 *
436  DO 200 i = 1, n
437  IF( rwork( i ).GT.safe2 ) THEN
438  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
439  ELSE
440  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
441  $ safe1
442  END IF
443  200 CONTINUE
444 *
445  kase = 0
446  210 CONTINUE
447  CALL clacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
448  IF( kase.NE.0 ) THEN
449  IF( kase.EQ.1 ) THEN
450 *
451 * Multiply by diag(W)*inv(op(A)**H).
452 *
453  CALL ctpsv( uplo, transt, diag, n, ap, work, 1 )
454  DO 220 i = 1, n
455  work( i ) = rwork( i )*work( i )
456  220 CONTINUE
457  ELSE
458 *
459 * Multiply by inv(op(A))*diag(W).
460 *
461  DO 230 i = 1, n
462  work( i ) = rwork( i )*work( i )
463  230 CONTINUE
464  CALL ctpsv( uplo, transn, diag, n, ap, work, 1 )
465  END IF
466  GO TO 210
467  END IF
468 *
469 * Normalize error.
470 *
471  lstres = zero
472  DO 240 i = 1, n
473  lstres = max( lstres, cabs1( x( i, j ) ) )
474  240 CONTINUE
475  IF( lstres.NE.zero )
476  $ ferr( j ) = ferr( j ) / lstres
477 *
478  250 CONTINUE
479 *
480  RETURN
481 *
482 * End of CTPRFS
483 *
ctpmv
subroutine ctpmv(UPLO, TRANS, DIAG, N, AP, X, INCX)
CTPMV
Definition: ctpmv.f:144
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
ctpsv
subroutine ctpsv(UPLO, TRANS, DIAG, N, AP, X, INCX)
CTPSV
Definition: ctpsv.f:146
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
caxpy
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:53
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
clacn2
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:135

Here is the call graph for this function:

Here is the caller graph for this function: