SRC/ztrrfs.f

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