SRC/ztbrfs.f

   1 *> \brief \b ZTBRFS
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZTBRFS + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztbrfs.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztbrfs.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztbrfs.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZTBRFS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
  22 *                          LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          DIAG, TRANS, UPLO
  26 *       INTEGER            INFO, KD, LDAB, LDB, LDX, N, NRHS
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * )
  30 *       COMPLEX*16         AB( LDAB, * ), B( LDB, * ), WORK( * ),
  31 *      $                   X( LDX, * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> ZTBRFS provides error bounds and backward error estimates for the
  41 *> solution to a system of linear equations with a triangular band
  42 *> coefficient matrix.
  43 *>
  44 *> The solution matrix X must be computed by ZTBTRS or some other
  45 *> means before entering this routine.  ZTBRFS does not do iterative
  46 *> refinement because doing so cannot improve the backward error.
  47 *> \endverbatim
  48 *
  49 *  Arguments:
  50 *  ==========
  51 *
  52 *> \param[in] UPLO
  53 *> \verbatim
  54 *>          UPLO is CHARACTER*1
  55 *>          = 'U':  A is upper triangular;
  56 *>          = 'L':  A is lower triangular.
  57 *> \endverbatim
  58 *>
  59 *> \param[in] TRANS
  60 *> \verbatim
  61 *>          TRANS is CHARACTER*1
  62 *>          Specifies the form of the system of equations:
  63 *>          = 'N':  A * X = B     (No transpose)
  64 *>          = 'T':  A**T * X = B  (Transpose)
  65 *>          = 'C':  A**H * X = B  (Conjugate transpose)
  66 *> \endverbatim
  67 *>
  68 *> \param[in] DIAG
  69 *> \verbatim
  70 *>          DIAG is CHARACTER*1
  71 *>          = 'N':  A is non-unit triangular;
  72 *>          = 'U':  A is unit triangular.
  73 *> \endverbatim
  74 *>
  75 *> \param[in] N
  76 *> \verbatim
  77 *>          N is INTEGER
  78 *>          The order of the matrix A.  N >= 0.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] KD
  82 *> \verbatim
  83 *>          KD is INTEGER
  84 *>          The number of superdiagonals or subdiagonals of the
  85 *>          triangular band matrix A.  KD >= 0.
  86 *> \endverbatim
  87 *>
  88 *> \param[in] NRHS
  89 *> \verbatim
  90 *>          NRHS is INTEGER
  91 *>          The number of right hand sides, i.e., the number of columns
  92 *>          of the matrices B and X.  NRHS >= 0.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] AB
  96 *> \verbatim
  97 *>          AB is COMPLEX*16 array, dimension (LDAB,N)
  98 *>          The upper or lower triangular band matrix A, stored in the
  99 *>          first kd+1 rows of the array. The j-th column of A is stored
 100 *>          in the j-th column of the array AB as follows:
 101 *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
 102 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
 103 *>          If DIAG = 'U', the diagonal elements of A are not referenced
 104 *>          and are assumed to be 1.
 105 *> \endverbatim
 106 *>
 107 *> \param[in] LDAB
 108 *> \verbatim
 109 *>          LDAB is INTEGER
 110 *>          The leading dimension of the array AB.  LDAB >= KD+1.
 111 *> \endverbatim
 112 *>
 113 *> \param[in] B
 114 *> \verbatim
 115 *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
 116 *>          The right hand side matrix B.
 117 *> \endverbatim
 118 *>
 119 *> \param[in] LDB
 120 *> \verbatim
 121 *>          LDB is INTEGER
 122 *>          The leading dimension of the array B.  LDB >= max(1,N).
 123 *> \endverbatim
 124 *>
 125 *> \param[in] X
 126 *> \verbatim
 127 *>          X is COMPLEX*16 array, dimension (LDX,NRHS)
 128 *>          The solution matrix X.
 129 *> \endverbatim
 130 *>
 131 *> \param[in] LDX
 132 *> \verbatim
 133 *>          LDX is INTEGER
 134 *>          The leading dimension of the array X.  LDX >= max(1,N).
 135 *> \endverbatim
 136 *>
 137 *> \param[out] FERR
 138 *> \verbatim
 139 *>          FERR is DOUBLE PRECISION array, dimension (NRHS)
 140 *>          The estimated forward error bound for each solution vector
 141 *>          X(j) (the j-th column of the solution matrix X).
 142 *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
 143 *>          is an estimated upper bound for the magnitude of the largest
 144 *>          element in (X(j) - XTRUE) divided by the magnitude of the
 145 *>          largest element in X(j).  The estimate is as reliable as
 146 *>          the estimate for RCOND, and is almost always a slight
 147 *>          overestimate of the true error.
 148 *> \endverbatim
 149 *>
 150 *> \param[out] BERR
 151 *> \verbatim
 152 *>          BERR is DOUBLE PRECISION array, dimension (NRHS)
 153 *>          The componentwise relative backward error of each solution
 154 *>          vector X(j) (i.e., the smallest relative change in
 155 *>          any element of A or B that makes X(j) an exact solution).
 156 *> \endverbatim
 157 *>
 158 *> \param[out] WORK
 159 *> \verbatim
 160 *>          WORK is COMPLEX*16 array, dimension (2*N)
 161 *> \endverbatim
 162 *>
 163 *> \param[out] RWORK
 164 *> \verbatim
 165 *>          RWORK is DOUBLE PRECISION array, dimension (N)
 166 *> \endverbatim
 167 *>
 168 *> \param[out] INFO
 169 *> \verbatim
 170 *>          INFO is INTEGER
 171 *>          = 0:  successful exit
 172 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 173 *> \endverbatim
 174 *
 175 *  Authors:
 176 *  ========
 177 *
 178 *> \author Univ. of Tennessee
 179 *> \author Univ. of California Berkeley
 180 *> \author Univ. of Colorado Denver
 181 *> \author NAG Ltd.
 182 *
 183 *> \date November 2011
 184 *
 185 *> \ingroup complex16OTHERcomputational
 186 *
 187 *  =====================================================================
 188       SUBROUTINE ZTBRFS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
 189      $                   LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
 190 *
 191 *  -- LAPACK computational routine (version 3.4.0) --
 192 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 193 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 194 *     November 2011
 195 *
 196 *     .. Scalar Arguments ..
 197       CHARACTER          DIAG, TRANS, UPLO
 198       INTEGER            INFO, KD, LDAB, LDB, LDX, N, NRHS
 199 *     ..
 200 *     .. Array Arguments ..
 201       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * )
 202       COMPLEX*16         AB( LDAB, * ), B( LDB, * ), WORK( * ),
 203      $                   X( LDX, * )
 204 *     ..
 205 *
 206 *  =====================================================================
 207 *
 208 *     .. Parameters ..
 209       DOUBLE PRECISION   ZERO
 210       PARAMETER          ( ZERO = 0.0D+0 )
 211       COMPLEX*16         ONE
 212       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
 213 *     ..
 214 *     .. Local Scalars ..
 215       LOGICAL            NOTRAN, NOUNIT, UPPER
 216       CHARACTER          TRANSN, TRANST
 217       INTEGER            I, J, K, KASE, NZ
 218       DOUBLE PRECISION   EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
 219       COMPLEX*16         ZDUM
 220 *     ..
 221 *     .. Local Arrays ..
 222       INTEGER            ISAVE( 3 )
 223 *     ..
 224 *     .. External Subroutines ..
 225       EXTERNAL           XERBLA, ZAXPY, ZCOPY, ZLACN2, ZTBMV, ZTBSV
 226 *     ..
 227 *     .. Intrinsic Functions ..
 228       INTRINSIC          ABS, DBLE, DIMAG, MAX, MIN
 229 *     ..
 230 *     .. External Functions ..
 231       LOGICAL            LSAME
 232       DOUBLE PRECISION   DLAMCH
 233       EXTERNAL           LSAME, DLAMCH
 234 *     ..
 235 *     .. Statement Functions ..
 236       DOUBLE PRECISION   CABS1
 237 *     ..
 238 *     .. Statement Function definitions ..
 239       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
 240 *     ..
 241 *     .. Executable Statements ..
 242 *
 243 *     Test the input parameters.
 244 *
 245       INFO = 0
 246       UPPER = LSAME( UPLO, 'U' )
 247       NOTRAN = LSAME( TRANS, 'N' )
 248       NOUNIT = LSAME( DIAG, 'N' )
 249 *
 250       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 251          INFO = -1
 252       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
 253      $         LSAME( TRANS, 'C' ) ) THEN
 254          INFO = -2
 255       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
 256          INFO = -3
 257       ELSE IF( N.LT.0 ) THEN
 258          INFO = -4
 259       ELSE IF( KD.LT.0 ) THEN
 260          INFO = -5
 261       ELSE IF( NRHS.LT.0 ) THEN
 262          INFO = -6
 263       ELSE IF( LDAB.LT.KD+1 ) THEN
 264          INFO = -8
 265       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 266          INFO = -10
 267       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
 268          INFO = -12
 269       END IF
 270       IF( INFO.NE.0 ) THEN
 271          CALL XERBLA( 'ZTBRFS', -INFO )
 272          RETURN
 273       END IF
 274 *
 275 *     Quick return if possible
 276 *
 277       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
 278          DO 10 J = 1, NRHS
 279             FERR( J ) = ZERO
 280             BERR( J ) = ZERO
 281    10    CONTINUE
 282          RETURN
 283       END IF
 284 *
 285       IF( NOTRAN ) THEN
 286          TRANSN = 'N'
 287          TRANST = 'C'
 288       ELSE
 289          TRANSN = 'C'
 290          TRANST = 'N'
 291       END IF
 292 *
 293 *     NZ = maximum number of nonzero elements in each row of A, plus 1
 294 *
 295       NZ = KD + 2
 296       EPS = DLAMCH( 'Epsilon' )
 297       SAFMIN = DLAMCH( 'Safe minimum' )
 298       SAFE1 = NZ*SAFMIN
 299       SAFE2 = SAFE1 / EPS
 300 *
 301 *     Do for each right hand side
 302 *
 303       DO 250 J = 1, NRHS
 304 *
 305 *        Compute residual R = B - op(A) * X,
 306 *        where op(A) = A, A**T, or A**H, depending on TRANS.
 307 *
 308          CALL ZCOPY( N, X( 1, J ), 1, WORK, 1 )
 309          CALL ZTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK, 1 )
 310          CALL ZAXPY( N, -ONE, B( 1, J ), 1, WORK, 1 )
 311 *
 312 *        Compute componentwise relative backward error from formula
 313 *
 314 *        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
 315 *
 316 *        where abs(Z) is the componentwise absolute value of the matrix
 317 *        or vector Z.  If the i-th component of the denominator is less
 318 *        than SAFE2, then SAFE1 is added to the i-th components of the
 319 *        numerator and denominator before dividing.
 320 *
 321          DO 20 I = 1, N
 322             RWORK( I ) = CABS1( B( I, J ) )
 323    20    CONTINUE
 324 *
 325          IF( NOTRAN ) THEN
 326 *
 327 *           Compute abs(A)*abs(X) + abs(B).
 328 *
 329             IF( UPPER ) THEN
 330                IF( NOUNIT ) THEN
 331                   DO 40 K = 1, N
 332                      XK = CABS1( X( K, J ) )
 333                      DO 30 I = MAX( 1, K-KD ), K
 334                         RWORK( I ) = RWORK( I ) +
 335      $                               CABS1( AB( KD+1+I-K, K ) )*XK
 336    30                CONTINUE
 337    40             CONTINUE
 338                ELSE
 339                   DO 60 K = 1, N
 340                      XK = CABS1( X( K, J ) )
 341                      DO 50 I = MAX( 1, K-KD ), K - 1
 342                         RWORK( I ) = RWORK( I ) +
 343      $                               CABS1( AB( KD+1+I-K, K ) )*XK
 344    50                CONTINUE
 345                      RWORK( K ) = RWORK( K ) + XK
 346    60             CONTINUE
 347                END IF
 348             ELSE
 349                IF( NOUNIT ) THEN
 350                   DO 80 K = 1, N
 351                      XK = CABS1( X( K, J ) )
 352                      DO 70 I = K, MIN( N, K+KD )
 353                         RWORK( I ) = RWORK( I ) +
 354      $                               CABS1( AB( 1+I-K, K ) )*XK
 355    70                CONTINUE
 356    80             CONTINUE
 357                ELSE
 358                   DO 100 K = 1, N
 359                      XK = CABS1( X( K, J ) )
 360                      DO 90 I = K + 1, MIN( N, K+KD )
 361                         RWORK( I ) = RWORK( I ) +
 362      $                               CABS1( AB( 1+I-K, K ) )*XK
 363    90                CONTINUE
 364                      RWORK( K ) = RWORK( K ) + XK
 365   100             CONTINUE
 366                END IF
 367             END IF
 368          ELSE
 369 *
 370 *           Compute abs(A**H)*abs(X) + abs(B).
 371 *
 372             IF( UPPER ) THEN
 373                IF( NOUNIT ) THEN
 374                   DO 120 K = 1, N
 375                      S = ZERO
 376                      DO 110 I = MAX( 1, K-KD ), K
 377                         S = S + CABS1( AB( KD+1+I-K, K ) )*
 378      $                      CABS1( X( I, J ) )
 379   110                CONTINUE
 380                      RWORK( K ) = RWORK( K ) + S
 381   120             CONTINUE
 382                ELSE
 383                   DO 140 K = 1, N
 384                      S = CABS1( X( K, J ) )
 385                      DO 130 I = MAX( 1, K-KD ), K - 1
 386                         S = S + CABS1( AB( KD+1+I-K, K ) )*
 387      $                      CABS1( X( I, J ) )
 388   130                CONTINUE
 389                      RWORK( K ) = RWORK( K ) + S
 390   140             CONTINUE
 391                END IF
 392             ELSE
 393                IF( NOUNIT ) THEN
 394                   DO 160 K = 1, N
 395                      S = ZERO
 396                      DO 150 I = K, MIN( N, K+KD )
 397                         S = S + CABS1( AB( 1+I-K, K ) )*
 398      $                      CABS1( X( I, J ) )
 399   150                CONTINUE
 400                      RWORK( K ) = RWORK( K ) + S
 401   160             CONTINUE
 402                ELSE
 403                   DO 180 K = 1, N
 404                      S = CABS1( X( K, J ) )
 405                      DO 170 I = K + 1, MIN( N, K+KD )
 406                         S = S + CABS1( AB( 1+I-K, K ) )*
 407      $                      CABS1( X( I, J ) )
 408   170                CONTINUE
 409                      RWORK( K ) = RWORK( K ) + S
 410   180             CONTINUE
 411                END IF
 412             END IF
 413          END IF
 414          S = ZERO
 415          DO 190 I = 1, N
 416             IF( RWORK( I ).GT.SAFE2 ) THEN
 417                S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
 418             ELSE
 419                S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
 420      $             ( RWORK( I )+SAFE1 ) )
 421             END IF
 422   190    CONTINUE
 423          BERR( J ) = S
 424 *
 425 *        Bound error from formula
 426 *
 427 *        norm(X - XTRUE) / norm(X) .le. FERR =
 428 *        norm( abs(inv(op(A)))*
 429 *           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
 430 *
 431 *        where
 432 *          norm(Z) is the magnitude of the largest component of Z
 433 *          inv(op(A)) is the inverse of op(A)
 434 *          abs(Z) is the componentwise absolute value of the matrix or
 435 *             vector Z
 436 *          NZ is the maximum number of nonzeros in any row of A, plus 1
 437 *          EPS is machine epsilon
 438 *
 439 *        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
 440 *        is incremented by SAFE1 if the i-th component of
 441 *        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
 442 *
 443 *        Use ZLACN2 to estimate the infinity-norm of the matrix
 444 *           inv(op(A)) * diag(W),
 445 *        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
 446 *
 447          DO 200 I = 1, N
 448             IF( RWORK( I ).GT.SAFE2 ) THEN
 449                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
 450             ELSE
 451                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
 452      $                      SAFE1
 453             END IF
 454   200    CONTINUE
 455 *
 456          KASE = 0
 457   210    CONTINUE
 458          CALL ZLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE )
 459          IF( KASE.NE.0 ) THEN
 460             IF( KASE.EQ.1 ) THEN
 461 *
 462 *              Multiply by diag(W)*inv(op(A)**H).
 463 *
 464                CALL ZTBSV( UPLO, TRANST, DIAG, N, KD, AB, LDAB, WORK,
 465      $                     1 )
 466                DO 220 I = 1, N
 467                   WORK( I ) = RWORK( I )*WORK( I )
 468   220          CONTINUE
 469             ELSE
 470 *
 471 *              Multiply by inv(op(A))*diag(W).
 472 *
 473                DO 230 I = 1, N
 474                   WORK( I ) = RWORK( I )*WORK( I )
 475   230          CONTINUE
 476                CALL ZTBSV( UPLO, TRANSN, DIAG, N, KD, AB, LDAB, WORK,
 477      $                     1 )
 478             END IF
 479             GO TO 210
 480          END IF
 481 *
 482 *        Normalize error.
 483 *
 484          LSTRES = ZERO
 485          DO 240 I = 1, N
 486             LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) )
 487   240    CONTINUE
 488          IF( LSTRES.NE.ZERO )
 489      $      FERR( J ) = FERR( J ) / LSTRES
 490 *
 491   250 CONTINUE
 492 *
 493       RETURN
 494 *
 495 *     End of ZTBRFS
 496 *
 497       END