SRC/ctprfs.f

   1 *> \brief \b CTPRFS
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CTPRFS + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctprfs.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctprfs.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctprfs.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CTPRFS( UPLO, TRANS, DIAG, N, NRHS, AP, B, LDB, X, LDX,
  22 *                          FERR, BERR, WORK, RWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          DIAG, TRANS, UPLO
  26 *       INTEGER            INFO, LDB, LDX, N, NRHS
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       REAL               BERR( * ), FERR( * ), RWORK( * )
  30 *       COMPLEX            AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> CTPRFS provides error bounds and backward error estimates for the
  40 *> solution to a system of linear equations with a triangular packed
  41 *> coefficient matrix.
  42 *>
  43 *> The solution matrix X must be computed by CTPTRS or some other
  44 *> means before entering this routine.  CTPRFS does not do iterative
  45 *> refinement because doing so cannot improve the backward error.
  46 *> \endverbatim
  47 *
  48 *  Arguments:
  49 *  ==========
  50 *
  51 *> \param[in] UPLO
  52 *> \verbatim
  53 *>          UPLO is CHARACTER*1
  54 *>          = 'U':  A is upper triangular;
  55 *>          = 'L':  A is lower triangular.
  56 *> \endverbatim
  57 *>
  58 *> \param[in] TRANS
  59 *> \verbatim
  60 *>          TRANS is CHARACTER*1
  61 *>          Specifies the form of the system of equations:
  62 *>          = 'N':  A * X = B     (No transpose)
  63 *>          = 'T':  A**T * X = B  (Transpose)
  64 *>          = 'C':  A**H * X = B  (Conjugate transpose)
  65 *> \endverbatim
  66 *>
  67 *> \param[in] DIAG
  68 *> \verbatim
  69 *>          DIAG is CHARACTER*1
  70 *>          = 'N':  A is non-unit triangular;
  71 *>          = 'U':  A is unit triangular.
  72 *> \endverbatim
  73 *>
  74 *> \param[in] N
  75 *> \verbatim
  76 *>          N is INTEGER
  77 *>          The order of the matrix A.  N >= 0.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] NRHS
  81 *> \verbatim
  82 *>          NRHS is INTEGER
  83 *>          The number of right hand sides, i.e., the number of columns
  84 *>          of the matrices B and X.  NRHS >= 0.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] AP
  88 *> \verbatim
  89 *>          AP is COMPLEX array, dimension (N*(N+1)/2)
  90 *>          The upper or lower triangular matrix A, packed columnwise in
  91 *>          a linear array.  The j-th column of A is stored in the array
  92 *>          AP as follows:
  93 *>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  94 *>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
  95 *>          If DIAG = 'U', the diagonal elements of A are not referenced
  96 *>          and are assumed to be 1.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] B
 100 *> \verbatim
 101 *>          B is COMPLEX array, dimension (LDB,NRHS)
 102 *>          The right hand side matrix B.
 103 *> \endverbatim
 104 *>
 105 *> \param[in] LDB
 106 *> \verbatim
 107 *>          LDB is INTEGER
 108 *>          The leading dimension of the array B.  LDB >= max(1,N).
 109 *> \endverbatim
 110 *>
 111 *> \param[in] X
 112 *> \verbatim
 113 *>          X is COMPLEX array, dimension (LDX,NRHS)
 114 *>          The solution matrix X.
 115 *> \endverbatim
 116 *>
 117 *> \param[in] LDX
 118 *> \verbatim
 119 *>          LDX is INTEGER
 120 *>          The leading dimension of the array X.  LDX >= max(1,N).
 121 *> \endverbatim
 122 *>
 123 *> \param[out] FERR
 124 *> \verbatim
 125 *>          FERR is REAL array, dimension (NRHS)
 126 *>          The estimated forward error bound for each solution vector
 127 *>          X(j) (the j-th column of the solution matrix X).
 128 *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
 129 *>          is an estimated upper bound for the magnitude of the largest
 130 *>          element in (X(j) - XTRUE) divided by the magnitude of the
 131 *>          largest element in X(j).  The estimate is as reliable as
 132 *>          the estimate for RCOND, and is almost always a slight
 133 *>          overestimate of the true error.
 134 *> \endverbatim
 135 *>
 136 *> \param[out] BERR
 137 *> \verbatim
 138 *>          BERR is REAL array, dimension (NRHS)
 139 *>          The componentwise relative backward error of each solution
 140 *>          vector X(j) (i.e., the smallest relative change in
 141 *>          any element of A or B that makes X(j) an exact solution).
 142 *> \endverbatim
 143 *>
 144 *> \param[out] WORK
 145 *> \verbatim
 146 *>          WORK is COMPLEX array, dimension (2*N)
 147 *> \endverbatim
 148 *>
 149 *> \param[out] RWORK
 150 *> \verbatim
 151 *>          RWORK is REAL array, dimension (N)
 152 *> \endverbatim
 153 *>
 154 *> \param[out] INFO
 155 *> \verbatim
 156 *>          INFO is INTEGER
 157 *>          = 0:  successful exit
 158 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 159 *> \endverbatim
 160 *
 161 *  Authors:
 162 *  ========
 163 *
 164 *> \author Univ. of Tennessee
 165 *> \author Univ. of California Berkeley
 166 *> \author Univ. of Colorado Denver
 167 *> \author NAG Ltd.
 168 *
 169 *> \date November 2011
 170 *
 171 *> \ingroup complexOTHERcomputational
 172 *
 173 *  =====================================================================
 174       SUBROUTINE CTPRFS( UPLO, TRANS, DIAG, N, NRHS, AP, B, LDB, X, LDX,
 175      $                   FERR, BERR, WORK, RWORK, INFO )
 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 *
 484       END