LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ ctprfs()

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.

Definition at line 172 of file ctprfs.f.

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