LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dtbrfs.f
Go to the documentation of this file.
00001 *> \brief \b DTBRFS
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DTBRFS + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtbrfs.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtbrfs.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtbrfs.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DTBRFS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
00022 *                          LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          DIAG, TRANS, UPLO
00026 *       INTEGER            INFO, KD, LDAB, LDB, LDX, N, NRHS
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            IWORK( * )
00030 *       DOUBLE PRECISION   AB( LDAB, * ), B( LDB, * ), BERR( * ),
00031 *      $                   FERR( * ), WORK( * ), X( LDX, * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> DTBRFS 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 DTBTRS or some other
00045 *> means before entering this routine.  DTBRFS 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 = 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (3*N)
00161 *> \endverbatim
00162 *>
00163 *> \param[out] IWORK
00164 *> \verbatim
00165 *>          IWORK is INTEGER 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 doubleOTHERcomputational
00186 *
00187 *  =====================================================================
00188       SUBROUTINE DTBRFS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
00189      $                   LDB, X, LDX, FERR, BERR, WORK, IWORK, 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       INTEGER            IWORK( * )
00202       DOUBLE PRECISION   AB( LDAB, * ), B( LDB, * ), BERR( * ),
00203      $                   FERR( * ), WORK( * ), X( LDX, * )
00204 *     ..
00205 *
00206 *  =====================================================================
00207 *
00208 *     .. Parameters ..
00209       DOUBLE PRECISION   ZERO
00210       PARAMETER          ( ZERO = 0.0D+0 )
00211       DOUBLE PRECISION   ONE
00212       PARAMETER          ( ONE = 1.0D+0 )
00213 *     ..
00214 *     .. Local Scalars ..
00215       LOGICAL            NOTRAN, NOUNIT, UPPER
00216       CHARACTER          TRANST
00217       INTEGER            I, J, K, KASE, NZ
00218       DOUBLE PRECISION   EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
00219 *     ..
00220 *     .. Local Arrays ..
00221       INTEGER            ISAVE( 3 )
00222 *     ..
00223 *     .. External Subroutines ..
00224       EXTERNAL           DAXPY, DCOPY, DLACN2, DTBMV, DTBSV, XERBLA
00225 *     ..
00226 *     .. Intrinsic Functions ..
00227       INTRINSIC          ABS, MAX, MIN
00228 *     ..
00229 *     .. External Functions ..
00230       LOGICAL            LSAME
00231       DOUBLE PRECISION   DLAMCH
00232       EXTERNAL           LSAME, DLAMCH
00233 *     ..
00234 *     .. Executable Statements ..
00235 *
00236 *     Test the input parameters.
00237 *
00238       INFO = 0
00239       UPPER = LSAME( UPLO, 'U' )
00240       NOTRAN = LSAME( TRANS, 'N' )
00241       NOUNIT = LSAME( DIAG, 'N' )
00242 *
00243       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00244          INFO = -1
00245       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
00246      $         LSAME( TRANS, 'C' ) ) THEN
00247          INFO = -2
00248       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
00249          INFO = -3
00250       ELSE IF( N.LT.0 ) THEN
00251          INFO = -4
00252       ELSE IF( KD.LT.0 ) THEN
00253          INFO = -5
00254       ELSE IF( NRHS.LT.0 ) THEN
00255          INFO = -6
00256       ELSE IF( LDAB.LT.KD+1 ) THEN
00257          INFO = -8
00258       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00259          INFO = -10
00260       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00261          INFO = -12
00262       END IF
00263       IF( INFO.NE.0 ) THEN
00264          CALL XERBLA( 'DTBRFS', -INFO )
00265          RETURN
00266       END IF
00267 *
00268 *     Quick return if possible
00269 *
00270       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
00271          DO 10 J = 1, NRHS
00272             FERR( J ) = ZERO
00273             BERR( J ) = ZERO
00274    10    CONTINUE
00275          RETURN
00276       END IF
00277 *
00278       IF( NOTRAN ) THEN
00279          TRANST = 'T'
00280       ELSE
00281          TRANST = 'N'
00282       END IF
00283 *
00284 *     NZ = maximum number of nonzero elements in each row of A, plus 1
00285 *
00286       NZ = KD + 2
00287       EPS = DLAMCH( 'Epsilon' )
00288       SAFMIN = DLAMCH( 'Safe minimum' )
00289       SAFE1 = NZ*SAFMIN
00290       SAFE2 = SAFE1 / EPS
00291 *
00292 *     Do for each right hand side
00293 *
00294       DO 250 J = 1, NRHS
00295 *
00296 *        Compute residual R = B - op(A) * X,
00297 *        where op(A) = A or A**T, depending on TRANS.
00298 *
00299          CALL DCOPY( N, X( 1, J ), 1, WORK( N+1 ), 1 )
00300          CALL DTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK( N+1 ),
00301      $               1 )
00302          CALL DAXPY( N, -ONE, B( 1, J ), 1, WORK( N+1 ), 1 )
00303 *
00304 *        Compute componentwise relative backward error from formula
00305 *
00306 *        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
00307 *
00308 *        where abs(Z) is the componentwise absolute value of the matrix
00309 *        or vector Z.  If the i-th component of the denominator is less
00310 *        than SAFE2, then SAFE1 is added to the i-th components of the
00311 *        numerator and denominator before dividing.
00312 *
00313          DO 20 I = 1, N
00314             WORK( I ) = ABS( B( I, J ) )
00315    20    CONTINUE
00316 *
00317          IF( NOTRAN ) THEN
00318 *
00319 *           Compute abs(A)*abs(X) + abs(B).
00320 *
00321             IF( UPPER ) THEN
00322                IF( NOUNIT ) THEN
00323                   DO 40 K = 1, N
00324                      XK = ABS( X( K, J ) )
00325                      DO 30 I = MAX( 1, K-KD ), K
00326                         WORK( I ) = WORK( I ) +
00327      $                              ABS( AB( KD+1+I-K, K ) )*XK
00328    30                CONTINUE
00329    40             CONTINUE
00330                ELSE
00331                   DO 60 K = 1, N
00332                      XK = ABS( X( K, J ) )
00333                      DO 50 I = MAX( 1, K-KD ), K - 1
00334                         WORK( I ) = WORK( I ) +
00335      $                              ABS( AB( KD+1+I-K, K ) )*XK
00336    50                CONTINUE
00337                      WORK( K ) = WORK( K ) + XK
00338    60             CONTINUE
00339                END IF
00340             ELSE
00341                IF( NOUNIT ) THEN
00342                   DO 80 K = 1, N
00343                      XK = ABS( X( K, J ) )
00344                      DO 70 I = K, MIN( N, K+KD )
00345                         WORK( I ) = WORK( I ) + ABS( AB( 1+I-K, K ) )*XK
00346    70                CONTINUE
00347    80             CONTINUE
00348                ELSE
00349                   DO 100 K = 1, N
00350                      XK = ABS( X( K, J ) )
00351                      DO 90 I = K + 1, MIN( N, K+KD )
00352                         WORK( I ) = WORK( I ) + ABS( AB( 1+I-K, K ) )*XK
00353    90                CONTINUE
00354                      WORK( K ) = WORK( K ) + XK
00355   100             CONTINUE
00356                END IF
00357             END IF
00358          ELSE
00359 *
00360 *           Compute abs(A**T)*abs(X) + abs(B).
00361 *
00362             IF( UPPER ) THEN
00363                IF( NOUNIT ) THEN
00364                   DO 120 K = 1, N
00365                      S = ZERO
00366                      DO 110 I = MAX( 1, K-KD ), K
00367                         S = S + ABS( AB( KD+1+I-K, K ) )*
00368      $                      ABS( X( I, J ) )
00369   110                CONTINUE
00370                      WORK( K ) = WORK( K ) + S
00371   120             CONTINUE
00372                ELSE
00373                   DO 140 K = 1, N
00374                      S = ABS( X( K, J ) )
00375                      DO 130 I = MAX( 1, K-KD ), K - 1
00376                         S = S + ABS( AB( KD+1+I-K, K ) )*
00377      $                      ABS( X( I, J ) )
00378   130                CONTINUE
00379                      WORK( K ) = WORK( K ) + S
00380   140             CONTINUE
00381                END IF
00382             ELSE
00383                IF( NOUNIT ) THEN
00384                   DO 160 K = 1, N
00385                      S = ZERO
00386                      DO 150 I = K, MIN( N, K+KD )
00387                         S = S + ABS( AB( 1+I-K, K ) )*ABS( X( I, J ) )
00388   150                CONTINUE
00389                      WORK( K ) = WORK( K ) + S
00390   160             CONTINUE
00391                ELSE
00392                   DO 180 K = 1, N
00393                      S = ABS( X( K, J ) )
00394                      DO 170 I = K + 1, MIN( N, K+KD )
00395                         S = S + ABS( AB( 1+I-K, K ) )*ABS( X( I, J ) )
00396   170                CONTINUE
00397                      WORK( K ) = WORK( K ) + S
00398   180             CONTINUE
00399                END IF
00400             END IF
00401          END IF
00402          S = ZERO
00403          DO 190 I = 1, N
00404             IF( WORK( I ).GT.SAFE2 ) THEN
00405                S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
00406             ELSE
00407                S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
00408      $             ( WORK( I )+SAFE1 ) )
00409             END IF
00410   190    CONTINUE
00411          BERR( J ) = S
00412 *
00413 *        Bound error from formula
00414 *
00415 *        norm(X - XTRUE) / norm(X) .le. FERR =
00416 *        norm( abs(inv(op(A)))*
00417 *           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
00418 *
00419 *        where
00420 *          norm(Z) is the magnitude of the largest component of Z
00421 *          inv(op(A)) is the inverse of op(A)
00422 *          abs(Z) is the componentwise absolute value of the matrix or
00423 *             vector Z
00424 *          NZ is the maximum number of nonzeros in any row of A, plus 1
00425 *          EPS is machine epsilon
00426 *
00427 *        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
00428 *        is incremented by SAFE1 if the i-th component of
00429 *        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
00430 *
00431 *        Use DLACN2 to estimate the infinity-norm of the matrix
00432 *           inv(op(A)) * diag(W),
00433 *        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
00434 *
00435          DO 200 I = 1, N
00436             IF( WORK( I ).GT.SAFE2 ) THEN
00437                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
00438             ELSE
00439                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
00440             END IF
00441   200    CONTINUE
00442 *
00443          KASE = 0
00444   210    CONTINUE
00445          CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
00446      $                KASE, ISAVE )
00447          IF( KASE.NE.0 ) THEN
00448             IF( KASE.EQ.1 ) THEN
00449 *
00450 *              Multiply by diag(W)*inv(op(A)**T).
00451 *
00452                CALL DTBSV( UPLO, TRANST, DIAG, N, KD, AB, LDAB,
00453      $                     WORK( N+1 ), 1 )
00454                DO 220 I = 1, N
00455                   WORK( N+I ) = WORK( I )*WORK( N+I )
00456   220          CONTINUE
00457             ELSE
00458 *
00459 *              Multiply by inv(op(A))*diag(W).
00460 *
00461                DO 230 I = 1, N
00462                   WORK( N+I ) = WORK( I )*WORK( N+I )
00463   230          CONTINUE
00464                CALL DTBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB,
00465      $                     WORK( N+1 ), 1 )
00466             END IF
00467             GO TO 210
00468          END IF
00469 *
00470 *        Normalize error.
00471 *
00472          LSTRES = ZERO
00473          DO 240 I = 1, N
00474             LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
00475   240    CONTINUE
00476          IF( LSTRES.NE.ZERO )
00477      $      FERR( J ) = FERR( J ) / LSTRES
00478 *
00479   250 CONTINUE
00480 *
00481       RETURN
00482 *
00483 *     End of DTBRFS
00484 *
00485       END
 All Files Functions