LAPACK: SRC/dtbrfs.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DTBRFS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
00002      $                   LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          DIAG, TRANS, UPLO
00013       INTEGER            INFO, KD, LDAB, LDB, LDX, N, NRHS
00014 *     ..
00015 *     .. Array Arguments ..
00016       INTEGER            IWORK( * )
00017       DOUBLE PRECISION   AB( LDAB, * ), B( LDB, * ), BERR( * ),
00018      $                   FERR( * ), WORK( * ), X( LDX, * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  DTBRFS provides error bounds and backward error estimates for the
00025 *  solution to a system of linear equations with a triangular band
00026 *  coefficient matrix.
00027 *
00028 *  The solution matrix X must be computed by DTBTRS or some other
00029 *  means before entering this routine.  DTBRFS does not do iterative
00030 *  refinement because doing so cannot improve the backward error.
00031 *
00032 *  Arguments
00033 *  =========
00034 *
00035 *  UPLO    (input) CHARACTER*1
00036 *          = 'U':  A is upper triangular;
00037 *          = 'L':  A is lower triangular.
00038 *
00039 *  TRANS   (input) CHARACTER*1
00040 *          Specifies the form of the system of equations:
00041 *          = 'N':  A * X = B  (No transpose)
00042 *          = 'T':  A**T * X = B  (Transpose)
00043 *          = 'C':  A**H * X = B  (Conjugate transpose = Transpose)
00044 *
00045 *  DIAG    (input) CHARACTER*1
00046 *          = 'N':  A is non-unit triangular;
00047 *          = 'U':  A is unit triangular.
00048 *
00049 *  N       (input) INTEGER
00050 *          The order of the matrix A.  N >= 0.
00051 *
00052 *  KD      (input) INTEGER
00053 *          The number of superdiagonals or subdiagonals of the
00054 *          triangular band matrix A.  KD >= 0.
00055 *
00056 *  NRHS    (input) INTEGER
00057 *          The number of right hand sides, i.e., the number of columns
00058 *          of the matrices B and X.  NRHS >= 0.
00059 *
00060 *  AB      (input) DOUBLE PRECISION array, dimension (LDAB,N)
00061 *          The upper or lower triangular band matrix A, stored in the
00062 *          first kd+1 rows of the array. The j-th column of A is stored
00063 *          in the j-th column of the array AB as follows:
00064 *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
00065 *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
00066 *          If DIAG = 'U', the diagonal elements of A are not referenced
00067 *          and are assumed to be 1.
00068 *
00069 *  LDAB    (input) INTEGER
00070 *          The leading dimension of the array AB.  LDAB >= KD+1.
00071 *
00072 *  B       (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
00073 *          The right hand side matrix B.
00074 *
00075 *  LDB     (input) INTEGER
00076 *          The leading dimension of the array B.  LDB >= max(1,N).
00077 *
00078 *  X       (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
00079 *          The solution matrix X.
00080 *
00081 *  LDX     (input) INTEGER
00082 *          The leading dimension of the array X.  LDX >= max(1,N).
00083 *
00084 *  FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
00085 *          The estimated forward error bound for each solution vector
00086 *          X(j) (the j-th column of the solution matrix X).
00087 *          If XTRUE is the true solution corresponding to X(j), FERR(j)
00088 *          is an estimated upper bound for the magnitude of the largest
00089 *          element in (X(j) - XTRUE) divided by the magnitude of the
00090 *          largest element in X(j).  The estimate is as reliable as
00091 *          the estimate for RCOND, and is almost always a slight
00092 *          overestimate of the true error.
00093 *
00094 *  BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
00095 *          The componentwise relative backward error of each solution
00096 *          vector X(j) (i.e., the smallest relative change in
00097 *          any element of A or B that makes X(j) an exact solution).
00098 *
00099 *  WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)
00100 *
00101 *  IWORK   (workspace) INTEGER array, dimension (N)
00102 *
00103 *  INFO    (output) INTEGER
00104 *          = 0:  successful exit
00105 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00106 *
00107 *  =====================================================================
00108 *
00109 *     .. Parameters ..
00110       DOUBLE PRECISION   ZERO
00111       PARAMETER          ( ZERO = 0.0D+0 )
00112       DOUBLE PRECISION   ONE
00113       PARAMETER          ( ONE = 1.0D+0 )
00114 *     ..
00115 *     .. Local Scalars ..
00116       LOGICAL            NOTRAN, NOUNIT, UPPER
00117       CHARACTER          TRANST
00118       INTEGER            I, J, K, KASE, NZ
00119       DOUBLE PRECISION   EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
00120 *     ..
00121 *     .. Local Arrays ..
00122       INTEGER            ISAVE( 3 )
00123 *     ..
00124 *     .. External Subroutines ..
00125       EXTERNAL           DAXPY, DCOPY, DLACN2, DTBMV, DTBSV, XERBLA
00126 *     ..
00127 *     .. Intrinsic Functions ..
00128       INTRINSIC          ABS, MAX, MIN
00129 *     ..
00130 *     .. External Functions ..
00131       LOGICAL            LSAME
00132       DOUBLE PRECISION   DLAMCH
00133       EXTERNAL           LSAME, DLAMCH
00134 *     ..
00135 *     .. Executable Statements ..
00136 *
00137 *     Test the input parameters.
00138 *
00139       INFO = 0
00140       UPPER = LSAME( UPLO, 'U' )
00141       NOTRAN = LSAME( TRANS, 'N' )
00142       NOUNIT = LSAME( DIAG, 'N' )
00143 *
00144       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00145          INFO = -1
00146       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
00147      $         LSAME( TRANS, 'C' ) ) THEN
00148          INFO = -2
00149       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
00150          INFO = -3
00151       ELSE IF( N.LT.0 ) THEN
00152          INFO = -4
00153       ELSE IF( KD.LT.0 ) THEN
00154          INFO = -5
00155       ELSE IF( NRHS.LT.0 ) THEN
00156          INFO = -6
00157       ELSE IF( LDAB.LT.KD+1 ) THEN
00158          INFO = -8
00159       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00160          INFO = -10
00161       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00162          INFO = -12
00163       END IF
00164       IF( INFO.NE.0 ) THEN
00165          CALL XERBLA( 'DTBRFS', -INFO )
00166          RETURN
00167       END IF
00168 *
00169 *     Quick return if possible
00170 *
00171       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
00172          DO 10 J = 1, NRHS
00173             FERR( J ) = ZERO
00174             BERR( J ) = ZERO
00175    10    CONTINUE
00176          RETURN
00177       END IF
00178 *
00179       IF( NOTRAN ) THEN
00180          TRANST = 'T'
00181       ELSE
00182          TRANST = 'N'
00183       END IF
00184 *
00185 *     NZ = maximum number of nonzero elements in each row of A, plus 1
00186 *
00187       NZ = KD + 2
00188       EPS = DLAMCH( 'Epsilon' )
00189       SAFMIN = DLAMCH( 'Safe minimum' )
00190       SAFE1 = NZ*SAFMIN
00191       SAFE2 = SAFE1 / EPS
00192 *
00193 *     Do for each right hand side
00194 *
00195       DO 250 J = 1, NRHS
00196 *
00197 *        Compute residual R = B - op(A) * X,
00198 *        where op(A) = A or A', depending on TRANS.
00199 *
00200          CALL DCOPY( N, X( 1, J ), 1, WORK( N+1 ), 1 )
00201          CALL DTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK( N+1 ),
00202      $               1 )
00203          CALL DAXPY( N, -ONE, B( 1, J ), 1, WORK( N+1 ), 1 )
00204 *
00205 *        Compute componentwise relative backward error from formula
00206 *
00207 *        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
00208 *
00209 *        where abs(Z) is the componentwise absolute value of the matrix
00210 *        or vector Z.  If the i-th component of the denominator is less
00211 *        than SAFE2, then SAFE1 is added to the i-th components of the
00212 *        numerator and denominator before dividing.
00213 *
00214          DO 20 I = 1, N
00215             WORK( I ) = ABS( B( I, J ) )
00216    20    CONTINUE
00217 *
00218          IF( NOTRAN ) THEN
00219 *
00220 *           Compute abs(A)*abs(X) + abs(B).
00221 *
00222             IF( UPPER ) THEN
00223                IF( NOUNIT ) THEN
00224                   DO 40 K = 1, N
00225                      XK = ABS( X( K, J ) )
00226                      DO 30 I = MAX( 1, K-KD ), K
00227                         WORK( I ) = WORK( I ) +
00228      $                              ABS( AB( KD+1+I-K, K ) )*XK
00229    30                CONTINUE
00230    40             CONTINUE
00231                ELSE
00232                   DO 60 K = 1, N
00233                      XK = ABS( X( K, J ) )
00234                      DO 50 I = MAX( 1, K-KD ), K - 1
00235                         WORK( I ) = WORK( I ) +
00236      $                              ABS( AB( KD+1+I-K, K ) )*XK
00237    50                CONTINUE
00238                      WORK( K ) = WORK( K ) + XK
00239    60             CONTINUE
00240                END IF
00241             ELSE
00242                IF( NOUNIT ) THEN
00243                   DO 80 K = 1, N
00244                      XK = ABS( X( K, J ) )
00245                      DO 70 I = K, MIN( N, K+KD )
00246                         WORK( I ) = WORK( I ) + ABS( AB( 1+I-K, K ) )*XK
00247    70                CONTINUE
00248    80             CONTINUE
00249                ELSE
00250                   DO 100 K = 1, N
00251                      XK = ABS( X( K, J ) )
00252                      DO 90 I = K + 1, MIN( N, K+KD )
00253                         WORK( I ) = WORK( I ) + ABS( AB( 1+I-K, K ) )*XK
00254    90                CONTINUE
00255                      WORK( K ) = WORK( K ) + XK
00256   100             CONTINUE
00257                END IF
00258             END IF
00259          ELSE
00260 *
00261 *           Compute abs(A')*abs(X) + abs(B).
00262 *
00263             IF( UPPER ) THEN
00264                IF( NOUNIT ) THEN
00265                   DO 120 K = 1, N
00266                      S = ZERO
00267                      DO 110 I = MAX( 1, K-KD ), K
00268                         S = S + ABS( AB( KD+1+I-K, K ) )*
00269      $                      ABS( X( I, J ) )
00270   110                CONTINUE
00271                      WORK( K ) = WORK( K ) + S
00272   120             CONTINUE
00273                ELSE
00274                   DO 140 K = 1, N
00275                      S = ABS( X( K, J ) )
00276                      DO 130 I = MAX( 1, K-KD ), K - 1
00277                         S = S + ABS( AB( KD+1+I-K, K ) )*
00278      $                      ABS( X( I, J ) )
00279   130                CONTINUE
00280                      WORK( K ) = WORK( K ) + S
00281   140             CONTINUE
00282                END IF
00283             ELSE
00284                IF( NOUNIT ) THEN
00285                   DO 160 K = 1, N
00286                      S = ZERO
00287                      DO 150 I = K, MIN( N, K+KD )
00288                         S = S + ABS( AB( 1+I-K, K ) )*ABS( X( I, J ) )
00289   150                CONTINUE
00290                      WORK( K ) = WORK( K ) + S
00291   160             CONTINUE
00292                ELSE
00293                   DO 180 K = 1, N
00294                      S = ABS( X( K, J ) )
00295                      DO 170 I = K + 1, MIN( N, K+KD )
00296                         S = S + ABS( AB( 1+I-K, K ) )*ABS( X( I, J ) )
00297   170                CONTINUE
00298                      WORK( K ) = WORK( K ) + S
00299   180             CONTINUE
00300                END IF
00301             END IF
00302          END IF
00303          S = ZERO
00304          DO 190 I = 1, N
00305             IF( WORK( I ).GT.SAFE2 ) THEN
00306                S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
00307             ELSE
00308                S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
00309      $             ( WORK( I )+SAFE1 ) )
00310             END IF
00311   190    CONTINUE
00312          BERR( J ) = S
00313 *
00314 *        Bound error from formula
00315 *
00316 *        norm(X - XTRUE) / norm(X) .le. FERR =
00317 *        norm( abs(inv(op(A)))*
00318 *           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
00319 *
00320 *        where
00321 *          norm(Z) is the magnitude of the largest component of Z
00322 *          inv(op(A)) is the inverse of op(A)
00323 *          abs(Z) is the componentwise absolute value of the matrix or
00324 *             vector Z
00325 *          NZ is the maximum number of nonzeros in any row of A, plus 1
00326 *          EPS is machine epsilon
00327 *
00328 *        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
00329 *        is incremented by SAFE1 if the i-th component of
00330 *        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
00331 *
00332 *        Use DLACN2 to estimate the infinity-norm of the matrix
00333 *           inv(op(A)) * diag(W),
00334 *        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
00335 *
00336          DO 200 I = 1, N
00337             IF( WORK( I ).GT.SAFE2 ) THEN
00338                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
00339             ELSE
00340                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
00341             END IF
00342   200    CONTINUE
00343 *
00344          KASE = 0
00345   210    CONTINUE
00346          CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
00347      $                KASE, ISAVE )
00348          IF( KASE.NE.0 ) THEN
00349             IF( KASE.EQ.1 ) THEN
00350 *
00351 *              Multiply by diag(W)*inv(op(A)').
00352 *
00353                CALL DTBSV( UPLO, TRANST, DIAG, N, KD, AB, LDAB,
00354      $                     WORK( N+1 ), 1 )
00355                DO 220 I = 1, N
00356                   WORK( N+I ) = WORK( I )*WORK( N+I )
00357   220          CONTINUE
00358             ELSE
00359 *
00360 *              Multiply by inv(op(A))*diag(W).
00361 *
00362                DO 230 I = 1, N
00363                   WORK( N+I ) = WORK( I )*WORK( N+I )
00364   230          CONTINUE
00365                CALL DTBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB,
00366      $                     WORK( N+1 ), 1 )
00367             END IF
00368             GO TO 210
00369          END IF
00370 *
00371 *        Normalize error.
00372 *
00373          LSTRES = ZERO
00374          DO 240 I = 1, N
00375             LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
00376   240    CONTINUE
00377          IF( LSTRES.NE.ZERO )
00378      $      FERR( J ) = FERR( J ) / LSTRES
00379 *
00380   250 CONTINUE
00381 *
00382       RETURN
00383 *
00384 *     End of DTBRFS
00385 *
00386       END