LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ctbrfs.f
Go to the documentation of this file.
00001 *> \brief \b CTBRFS
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CTBRFS + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctbrfs.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctbrfs.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctbrfs.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CTBRFS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
00022 *                          LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          DIAG, TRANS, UPLO
00026 *       INTEGER            INFO, KD, LDAB, LDB, LDX, N, NRHS
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       REAL               BERR( * ), FERR( * ), RWORK( * )
00030 *       COMPLEX            AB( LDAB, * ), B( LDB, * ), WORK( * ),
00031 *      $                   X( LDX, * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> CTBRFS provides error bounds and backward error estimates for the
00041 *> solution to a system of linear equations with a triangular band
00042 *> coefficient matrix.
00043 *>
00044 *> The solution matrix X must be computed by CTBTRS or some other
00045 *> means before entering this routine.  CTBRFS does not do iterative
00046 *> refinement because doing so cannot improve the backward error.
00047 *> \endverbatim
00048 *
00049 *  Arguments:
00050 *  ==========
00051 *
00052 *> \param[in] UPLO
00053 *> \verbatim
00054 *>          UPLO is CHARACTER*1
00055 *>          = 'U':  A is upper triangular;
00056 *>          = 'L':  A is lower triangular.
00057 *> \endverbatim
00058 *>
00059 *> \param[in] TRANS
00060 *> \verbatim
00061 *>          TRANS is CHARACTER*1
00062 *>          Specifies the form of the system of equations:
00063 *>          = 'N':  A * X = B     (No transpose)
00064 *>          = 'T':  A**T * X = B  (Transpose)
00065 *>          = 'C':  A**H * X = B  (Conjugate transpose)
00066 *> \endverbatim
00067 *>
00068 *> \param[in] DIAG
00069 *> \verbatim
00070 *>          DIAG is CHARACTER*1
00071 *>          = 'N':  A is non-unit triangular;
00072 *>          = 'U':  A is unit triangular.
00073 *> \endverbatim
00074 *>
00075 *> \param[in] N
00076 *> \verbatim
00077 *>          N is INTEGER
00078 *>          The order of the matrix A.  N >= 0.
00079 *> \endverbatim
00080 *>
00081 *> \param[in] KD
00082 *> \verbatim
00083 *>          KD is INTEGER
00084 *>          The number of superdiagonals or subdiagonals of the
00085 *>          triangular band matrix A.  KD >= 0.
00086 *> \endverbatim
00087 *>
00088 *> \param[in] NRHS
00089 *> \verbatim
00090 *>          NRHS is INTEGER
00091 *>          The number of right hand sides, i.e., the number of columns
00092 *>          of the matrices B and X.  NRHS >= 0.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] AB
00096 *> \verbatim
00097 *>          AB is COMPLEX array, dimension (LDAB,N)
00098 *>          The upper or lower triangular band matrix A, stored in the
00099 *>          first kd+1 rows of the array. The j-th column of A is stored
00100 *>          in the j-th column of the array AB as follows:
00101 *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
00102 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
00103 *>          If DIAG = 'U', the diagonal elements of A are not referenced
00104 *>          and are assumed to be 1.
00105 *> \endverbatim
00106 *>
00107 *> \param[in] LDAB
00108 *> \verbatim
00109 *>          LDAB is INTEGER
00110 *>          The leading dimension of the array AB.  LDAB >= KD+1.
00111 *> \endverbatim
00112 *>
00113 *> \param[in] B
00114 *> \verbatim
00115 *>          B is COMPLEX array, dimension (LDB,NRHS)
00116 *>          The right hand side matrix B.
00117 *> \endverbatim
00118 *>
00119 *> \param[in] LDB
00120 *> \verbatim
00121 *>          LDB is INTEGER
00122 *>          The leading dimension of the array B.  LDB >= max(1,N).
00123 *> \endverbatim
00124 *>
00125 *> \param[in] X
00126 *> \verbatim
00127 *>          X is COMPLEX array, dimension (LDX,NRHS)
00128 *>          The solution matrix X.
00129 *> \endverbatim
00130 *>
00131 *> \param[in] LDX
00132 *> \verbatim
00133 *>          LDX is INTEGER
00134 *>          The leading dimension of the array X.  LDX >= max(1,N).
00135 *> \endverbatim
00136 *>
00137 *> \param[out] FERR
00138 *> \verbatim
00139 *>          FERR is REAL array, dimension (NRHS)
00140 *>          The estimated forward error bound for each solution vector
00141 *>          X(j) (the j-th column of the solution matrix X).
00142 *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
00143 *>          is an estimated upper bound for the magnitude of the largest
00144 *>          element in (X(j) - XTRUE) divided by the magnitude of the
00145 *>          largest element in X(j).  The estimate is as reliable as
00146 *>          the estimate for RCOND, and is almost always a slight
00147 *>          overestimate of the true error.
00148 *> \endverbatim
00149 *>
00150 *> \param[out] BERR
00151 *> \verbatim
00152 *>          BERR is REAL array, dimension (NRHS)
00153 *>          The componentwise relative backward error of each solution
00154 *>          vector X(j) (i.e., the smallest relative change in
00155 *>          any element of A or B that makes X(j) an exact solution).
00156 *> \endverbatim
00157 *>
00158 *> \param[out] WORK
00159 *> \verbatim
00160 *>          WORK is COMPLEX array, dimension (2*N)
00161 *> \endverbatim
00162 *>
00163 *> \param[out] RWORK
00164 *> \verbatim
00165 *>          RWORK is REAL array, dimension (N)
00166 *> \endverbatim
00167 *>
00168 *> \param[out] INFO
00169 *> \verbatim
00170 *>          INFO is INTEGER
00171 *>          = 0:  successful exit
00172 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00173 *> \endverbatim
00174 *
00175 *  Authors:
00176 *  ========
00177 *
00178 *> \author Univ. of Tennessee 
00179 *> \author Univ. of California Berkeley 
00180 *> \author Univ. of Colorado Denver 
00181 *> \author NAG Ltd. 
00182 *
00183 *> \date November 2011
00184 *
00185 *> \ingroup complexOTHERcomputational
00186 *
00187 *  =====================================================================
00188       SUBROUTINE CTBRFS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
00189      $                   LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
00190 *
00191 *  -- LAPACK computational routine (version 3.4.0) --
00192 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00193 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00194 *     November 2011
00195 *
00196 *     .. Scalar Arguments ..
00197       CHARACTER          DIAG, TRANS, UPLO
00198       INTEGER            INFO, KD, LDAB, LDB, LDX, N, NRHS
00199 *     ..
00200 *     .. Array Arguments ..
00201       REAL               BERR( * ), FERR( * ), RWORK( * )
00202       COMPLEX            AB( LDAB, * ), B( LDB, * ), WORK( * ),
00203      $                   X( LDX, * )
00204 *     ..
00205 *
00206 *  =====================================================================
00207 *
00208 *     .. Parameters ..
00209       REAL               ZERO
00210       PARAMETER          ( ZERO = 0.0E+0 )
00211       COMPLEX            ONE
00212       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ) )
00213 *     ..
00214 *     .. Local Scalars ..
00215       LOGICAL            NOTRAN, NOUNIT, UPPER
00216       CHARACTER          TRANSN, TRANST
00217       INTEGER            I, J, K, KASE, NZ
00218       REAL               EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
00219       COMPLEX            ZDUM
00220 *     ..
00221 *     .. Local Arrays ..
00222       INTEGER            ISAVE( 3 )
00223 *     ..
00224 *     .. External Subroutines ..
00225       EXTERNAL           CAXPY, CCOPY, CLACN2, CTBMV, CTBSV, XERBLA
00226 *     ..
00227 *     .. Intrinsic Functions ..
00228       INTRINSIC          ABS, AIMAG, MAX, MIN, REAL
00229 *     ..
00230 *     .. External Functions ..
00231       LOGICAL            LSAME
00232       REAL               SLAMCH
00233       EXTERNAL           LSAME, SLAMCH
00234 *     ..
00235 *     .. Statement Functions ..
00236       REAL               CABS1
00237 *     ..
00238 *     .. Statement Function definitions ..
00239       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
00240 *     ..
00241 *     .. Executable Statements ..
00242 *
00243 *     Test the input parameters.
00244 *
00245       INFO = 0
00246       UPPER = LSAME( UPLO, 'U' )
00247       NOTRAN = LSAME( TRANS, 'N' )
00248       NOUNIT = LSAME( DIAG, 'N' )
00249 *
00250       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00251          INFO = -1
00252       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
00253      $         LSAME( TRANS, 'C' ) ) THEN
00254          INFO = -2
00255       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
00256          INFO = -3
00257       ELSE IF( N.LT.0 ) THEN
00258          INFO = -4
00259       ELSE IF( KD.LT.0 ) THEN
00260          INFO = -5
00261       ELSE IF( NRHS.LT.0 ) THEN
00262          INFO = -6
00263       ELSE IF( LDAB.LT.KD+1 ) THEN
00264          INFO = -8
00265       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00266          INFO = -10
00267       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00268          INFO = -12
00269       END IF
00270       IF( INFO.NE.0 ) THEN
00271          CALL XERBLA( 'CTBRFS', -INFO )
00272          RETURN
00273       END IF
00274 *
00275 *     Quick return if possible
00276 *
00277       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
00278          DO 10 J = 1, NRHS
00279             FERR( J ) = ZERO
00280             BERR( J ) = ZERO
00281    10    CONTINUE
00282          RETURN
00283       END IF
00284 *
00285       IF( NOTRAN ) THEN
00286          TRANSN = 'N'
00287          TRANST = 'C'
00288       ELSE
00289          TRANSN = 'C'
00290          TRANST = 'N'
00291       END IF
00292 *
00293 *     NZ = maximum number of nonzero elements in each row of A, plus 1
00294 *
00295       NZ = KD + 2
00296       EPS = SLAMCH( 'Epsilon' )
00297       SAFMIN = SLAMCH( 'Safe minimum' )
00298       SAFE1 = NZ*SAFMIN
00299       SAFE2 = SAFE1 / EPS
00300 *
00301 *     Do for each right hand side
00302 *
00303       DO 250 J = 1, NRHS
00304 *
00305 *        Compute residual R = B - op(A) * X,
00306 *        where op(A) = A, A**T, or A**H, depending on TRANS.
00307 *
00308          CALL CCOPY( N, X( 1, J ), 1, WORK, 1 )
00309          CALL CTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK, 1 )
00310          CALL CAXPY( N, -ONE, B( 1, J ), 1, WORK, 1 )
00311 *
00312 *        Compute componentwise relative backward error from formula
00313 *
00314 *        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
00315 *
00316 *        where abs(Z) is the componentwise absolute value of the matrix
00317 *        or vector Z.  If the i-th component of the denominator is less
00318 *        than SAFE2, then SAFE1 is added to the i-th components of the
00319 *        numerator and denominator before dividing.
00320 *
00321          DO 20 I = 1, N
00322             RWORK( I ) = CABS1( B( I, J ) )
00323    20    CONTINUE
00324 *
00325          IF( NOTRAN ) THEN
00326 *
00327 *           Compute abs(A)*abs(X) + abs(B).
00328 *
00329             IF( UPPER ) THEN
00330                IF( NOUNIT ) THEN
00331                   DO 40 K = 1, N
00332                      XK = CABS1( X( K, J ) )
00333                      DO 30 I = MAX( 1, K-KD ), K
00334                         RWORK( I ) = RWORK( I ) +
00335      $                               CABS1( AB( KD+1+I-K, K ) )*XK
00336    30                CONTINUE
00337    40             CONTINUE
00338                ELSE
00339                   DO 60 K = 1, N
00340                      XK = CABS1( X( K, J ) )
00341                      DO 50 I = MAX( 1, K-KD ), K - 1
00342                         RWORK( I ) = RWORK( I ) +
00343      $                               CABS1( AB( KD+1+I-K, K ) )*XK
00344    50                CONTINUE
00345                      RWORK( K ) = RWORK( K ) + XK
00346    60             CONTINUE
00347                END IF
00348             ELSE
00349                IF( NOUNIT ) THEN
00350                   DO 80 K = 1, N
00351                      XK = CABS1( X( K, J ) )
00352                      DO 70 I = K, MIN( N, K+KD )
00353                         RWORK( I ) = RWORK( I ) +
00354      $                               CABS1( AB( 1+I-K, K ) )*XK
00355    70                CONTINUE
00356    80             CONTINUE
00357                ELSE
00358                   DO 100 K = 1, N
00359                      XK = CABS1( X( K, J ) )
00360                      DO 90 I = K + 1, MIN( N, K+KD )
00361                         RWORK( I ) = RWORK( I ) +
00362      $                               CABS1( AB( 1+I-K, K ) )*XK
00363    90                CONTINUE
00364                      RWORK( K ) = RWORK( K ) + XK
00365   100             CONTINUE
00366                END IF
00367             END IF
00368          ELSE
00369 *
00370 *           Compute abs(A**H)*abs(X) + abs(B).
00371 *
00372             IF( UPPER ) THEN
00373                IF( NOUNIT ) THEN
00374                   DO 120 K = 1, N
00375                      S = ZERO
00376                      DO 110 I = MAX( 1, K-KD ), K
00377                         S = S + CABS1( AB( KD+1+I-K, K ) )*
00378      $                      CABS1( X( I, J ) )
00379   110                CONTINUE
00380                      RWORK( K ) = RWORK( K ) + S
00381   120             CONTINUE
00382                ELSE
00383                   DO 140 K = 1, N
00384                      S = CABS1( X( K, J ) )
00385                      DO 130 I = MAX( 1, K-KD ), K - 1
00386                         S = S + CABS1( AB( KD+1+I-K, K ) )*
00387      $                      CABS1( X( I, J ) )
00388   130                CONTINUE
00389                      RWORK( K ) = RWORK( K ) + S
00390   140             CONTINUE
00391                END IF
00392             ELSE
00393                IF( NOUNIT ) THEN
00394                   DO 160 K = 1, N
00395                      S = ZERO
00396                      DO 150 I = K, MIN( N, K+KD )
00397                         S = S + CABS1( AB( 1+I-K, K ) )*
00398      $                      CABS1( X( I, J ) )
00399   150                CONTINUE
00400                      RWORK( K ) = RWORK( K ) + S
00401   160             CONTINUE
00402                ELSE
00403                   DO 180 K = 1, N
00404                      S = CABS1( X( K, J ) )
00405                      DO 170 I = K + 1, MIN( N, K+KD )
00406                         S = S + CABS1( AB( 1+I-K, K ) )*
00407      $                      CABS1( X( I, J ) )
00408   170                CONTINUE
00409                      RWORK( K ) = RWORK( K ) + S
00410   180             CONTINUE
00411                END IF
00412             END IF
00413          END IF
00414          S = ZERO
00415          DO 190 I = 1, N
00416             IF( RWORK( I ).GT.SAFE2 ) THEN
00417                S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
00418             ELSE
00419                S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
00420      $             ( RWORK( I )+SAFE1 ) )
00421             END IF
00422   190    CONTINUE
00423          BERR( J ) = S
00424 *
00425 *        Bound error from formula
00426 *
00427 *        norm(X - XTRUE) / norm(X) .le. FERR =
00428 *        norm( abs(inv(op(A)))*
00429 *           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
00430 *
00431 *        where
00432 *          norm(Z) is the magnitude of the largest component of Z
00433 *          inv(op(A)) is the inverse of op(A)
00434 *          abs(Z) is the componentwise absolute value of the matrix or
00435 *             vector Z
00436 *          NZ is the maximum number of nonzeros in any row of A, plus 1
00437 *          EPS is machine epsilon
00438 *
00439 *        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
00440 *        is incremented by SAFE1 if the i-th component of
00441 *        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
00442 *
00443 *        Use CLACN2 to estimate the infinity-norm of the matrix
00444 *           inv(op(A)) * diag(W),
00445 *        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
00446 *
00447          DO 200 I = 1, N
00448             IF( RWORK( I ).GT.SAFE2 ) THEN
00449                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
00450             ELSE
00451                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
00452      $                      SAFE1
00453             END IF
00454   200    CONTINUE
00455 *
00456          KASE = 0
00457   210    CONTINUE
00458          CALL CLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE )
00459          IF( KASE.NE.0 ) THEN
00460             IF( KASE.EQ.1 ) THEN
00461 *
00462 *              Multiply by diag(W)*inv(op(A)**H).
00463 *
00464                CALL CTBSV( UPLO, TRANST, DIAG, N, KD, AB, LDAB, WORK,
00465      $                     1 )
00466                DO 220 I = 1, N
00467                   WORK( I ) = RWORK( I )*WORK( I )
00468   220          CONTINUE
00469             ELSE
00470 *
00471 *              Multiply by inv(op(A))*diag(W).
00472 *
00473                DO 230 I = 1, N
00474                   WORK( I ) = RWORK( I )*WORK( I )
00475   230          CONTINUE
00476                CALL CTBSV( UPLO, TRANSN, DIAG, N, KD, AB, LDAB, WORK,
00477      $                     1 )
00478             END IF
00479             GO TO 210
00480          END IF
00481 *
00482 *        Normalize error.
00483 *
00484          LSTRES = ZERO
00485          DO 240 I = 1, N
00486             LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) )
00487   240    CONTINUE
00488          IF( LSTRES.NE.ZERO )
00489      $      FERR( J ) = FERR( J ) / LSTRES
00490 *
00491   250 CONTINUE
00492 *
00493       RETURN
00494 *
00495 *     End of CTBRFS
00496 *
00497       END
 All Files Functions