LAPACK: SRC/ctbrfs.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CTBRFS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
00002      $                   LDB, X, LDX, FERR, BERR, WORK, RWORK, 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 CLACN2 in place of CLACON, 10 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       REAL               BERR( * ), FERR( * ), RWORK( * )
00017       COMPLEX            AB( LDAB, * ), B( LDB, * ), WORK( * ),
00018      $                   X( LDX, * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  CTBRFS 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 CTBTRS or some other
00029 *  means before entering this routine.  CTBRFS 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)
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) COMPLEX 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) COMPLEX 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) COMPLEX 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) REAL 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) REAL 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) COMPLEX array, dimension (2*N)
00100 *
00101 *  RWORK   (workspace) REAL 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       REAL               ZERO
00111       PARAMETER          ( ZERO = 0.0E+0 )
00112       COMPLEX            ONE
00113       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ) )
00114 *     ..
00115 *     .. Local Scalars ..
00116       LOGICAL            NOTRAN, NOUNIT, UPPER
00117       CHARACTER          TRANSN, TRANST
00118       INTEGER            I, J, K, KASE, NZ
00119       REAL               EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
00120       COMPLEX            ZDUM
00121 *     ..
00122 *     .. Local Arrays ..
00123       INTEGER            ISAVE( 3 )
00124 *     ..
00125 *     .. External Subroutines ..
00126       EXTERNAL           CAXPY, CCOPY, CLACN2, CTBMV, CTBSV, XERBLA
00127 *     ..
00128 *     .. Intrinsic Functions ..
00129       INTRINSIC          ABS, AIMAG, MAX, MIN, REAL
00130 *     ..
00131 *     .. External Functions ..
00132       LOGICAL            LSAME
00133       REAL               SLAMCH
00134       EXTERNAL           LSAME, SLAMCH
00135 *     ..
00136 *     .. Statement Functions ..
00137       REAL               CABS1
00138 *     ..
00139 *     .. Statement Function definitions ..
00140       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
00141 *     ..
00142 *     .. Executable Statements ..
00143 *
00144 *     Test the input parameters.
00145 *
00146       INFO = 0
00147       UPPER = LSAME( UPLO, 'U' )
00148       NOTRAN = LSAME( TRANS, 'N' )
00149       NOUNIT = LSAME( DIAG, 'N' )
00150 *
00151       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00152          INFO = -1
00153       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
00154      $         LSAME( TRANS, 'C' ) ) THEN
00155          INFO = -2
00156       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
00157          INFO = -3
00158       ELSE IF( N.LT.0 ) THEN
00159          INFO = -4
00160       ELSE IF( KD.LT.0 ) THEN
00161          INFO = -5
00162       ELSE IF( NRHS.LT.0 ) THEN
00163          INFO = -6
00164       ELSE IF( LDAB.LT.KD+1 ) THEN
00165          INFO = -8
00166       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00167          INFO = -10
00168       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00169          INFO = -12
00170       END IF
00171       IF( INFO.NE.0 ) THEN
00172          CALL XERBLA( 'CTBRFS', -INFO )
00173          RETURN
00174       END IF
00175 *
00176 *     Quick return if possible
00177 *
00178       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
00179          DO 10 J = 1, NRHS
00180             FERR( J ) = ZERO
00181             BERR( J ) = ZERO
00182    10    CONTINUE
00183          RETURN
00184       END IF
00185 *
00186       IF( NOTRAN ) THEN
00187          TRANSN = 'N'
00188          TRANST = 'C'
00189       ELSE
00190          TRANSN = 'C'
00191          TRANST = 'N'
00192       END IF
00193 *
00194 *     NZ = maximum number of nonzero elements in each row of A, plus 1
00195 *
00196       NZ = KD + 2
00197       EPS = SLAMCH( 'Epsilon' )
00198       SAFMIN = SLAMCH( 'Safe minimum' )
00199       SAFE1 = NZ*SAFMIN
00200       SAFE2 = SAFE1 / EPS
00201 *
00202 *     Do for each right hand side
00203 *
00204       DO 250 J = 1, NRHS
00205 *
00206 *        Compute residual R = B - op(A) * X,
00207 *        where op(A) = A, A**T, or A**H, depending on TRANS.
00208 *
00209          CALL CCOPY( N, X( 1, J ), 1, WORK, 1 )
00210          CALL CTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK, 1 )
00211          CALL CAXPY( N, -ONE, B( 1, J ), 1, WORK, 1 )
00212 *
00213 *        Compute componentwise relative backward error from formula
00214 *
00215 *        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
00216 *
00217 *        where abs(Z) is the componentwise absolute value of the matrix
00218 *        or vector Z.  If the i-th component of the denominator is less
00219 *        than SAFE2, then SAFE1 is added to the i-th components of the
00220 *        numerator and denominator before dividing.
00221 *
00222          DO 20 I = 1, N
00223             RWORK( I ) = CABS1( B( I, J ) )
00224    20    CONTINUE
00225 *
00226          IF( NOTRAN ) THEN
00227 *
00228 *           Compute abs(A)*abs(X) + abs(B).
00229 *
00230             IF( UPPER ) THEN
00231                IF( NOUNIT ) THEN
00232                   DO 40 K = 1, N
00233                      XK = CABS1( X( K, J ) )
00234                      DO 30 I = MAX( 1, K-KD ), K
00235                         RWORK( I ) = RWORK( I ) +
00236      $                               CABS1( AB( KD+1+I-K, K ) )*XK
00237    30                CONTINUE
00238    40             CONTINUE
00239                ELSE
00240                   DO 60 K = 1, N
00241                      XK = CABS1( X( K, J ) )
00242                      DO 50 I = MAX( 1, K-KD ), K - 1
00243                         RWORK( I ) = RWORK( I ) +
00244      $                               CABS1( AB( KD+1+I-K, K ) )*XK
00245    50                CONTINUE
00246                      RWORK( K ) = RWORK( K ) + XK
00247    60             CONTINUE
00248                END IF
00249             ELSE
00250                IF( NOUNIT ) THEN
00251                   DO 80 K = 1, N
00252                      XK = CABS1( X( K, J ) )
00253                      DO 70 I = K, MIN( N, K+KD )
00254                         RWORK( I ) = RWORK( I ) +
00255      $                               CABS1( AB( 1+I-K, K ) )*XK
00256    70                CONTINUE
00257    80             CONTINUE
00258                ELSE
00259                   DO 100 K = 1, N
00260                      XK = CABS1( X( K, J ) )
00261                      DO 90 I = K + 1, MIN( N, K+KD )
00262                         RWORK( I ) = RWORK( I ) +
00263      $                               CABS1( AB( 1+I-K, K ) )*XK
00264    90                CONTINUE
00265                      RWORK( K ) = RWORK( K ) + XK
00266   100             CONTINUE
00267                END IF
00268             END IF
00269          ELSE
00270 *
00271 *           Compute abs(A**H)*abs(X) + abs(B).
00272 *
00273             IF( UPPER ) THEN
00274                IF( NOUNIT ) THEN
00275                   DO 120 K = 1, N
00276                      S = ZERO
00277                      DO 110 I = MAX( 1, K-KD ), K
00278                         S = S + CABS1( AB( KD+1+I-K, K ) )*
00279      $                      CABS1( X( I, J ) )
00280   110                CONTINUE
00281                      RWORK( K ) = RWORK( K ) + S
00282   120             CONTINUE
00283                ELSE
00284                   DO 140 K = 1, N
00285                      S = CABS1( X( K, J ) )
00286                      DO 130 I = MAX( 1, K-KD ), K - 1
00287                         S = S + CABS1( AB( KD+1+I-K, K ) )*
00288      $                      CABS1( X( I, J ) )
00289   130                CONTINUE
00290                      RWORK( K ) = RWORK( K ) + S
00291   140             CONTINUE
00292                END IF
00293             ELSE
00294                IF( NOUNIT ) THEN
00295                   DO 160 K = 1, N
00296                      S = ZERO
00297                      DO 150 I = K, MIN( N, K+KD )
00298                         S = S + CABS1( AB( 1+I-K, K ) )*
00299      $                      CABS1( X( I, J ) )
00300   150                CONTINUE
00301                      RWORK( K ) = RWORK( K ) + S
00302   160             CONTINUE
00303                ELSE
00304                   DO 180 K = 1, N
00305                      S = CABS1( X( K, J ) )
00306                      DO 170 I = K + 1, MIN( N, K+KD )
00307                         S = S + CABS1( AB( 1+I-K, K ) )*
00308      $                      CABS1( X( I, J ) )
00309   170                CONTINUE
00310                      RWORK( K ) = RWORK( K ) + S
00311   180             CONTINUE
00312                END IF
00313             END IF
00314          END IF
00315          S = ZERO
00316          DO 190 I = 1, N
00317             IF( RWORK( I ).GT.SAFE2 ) THEN
00318                S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
00319             ELSE
00320                S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
00321      $             ( RWORK( I )+SAFE1 ) )
00322             END IF
00323   190    CONTINUE
00324          BERR( J ) = S
00325 *
00326 *        Bound error from formula
00327 *
00328 *        norm(X - XTRUE) / norm(X) .le. FERR =
00329 *        norm( abs(inv(op(A)))*
00330 *           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
00331 *
00332 *        where
00333 *          norm(Z) is the magnitude of the largest component of Z
00334 *          inv(op(A)) is the inverse of op(A)
00335 *          abs(Z) is the componentwise absolute value of the matrix or
00336 *             vector Z
00337 *          NZ is the maximum number of nonzeros in any row of A, plus 1
00338 *          EPS is machine epsilon
00339 *
00340 *        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
00341 *        is incremented by SAFE1 if the i-th component of
00342 *        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
00343 *
00344 *        Use CLACN2 to estimate the infinity-norm of the matrix
00345 *           inv(op(A)) * diag(W),
00346 *        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
00347 *
00348          DO 200 I = 1, N
00349             IF( RWORK( I ).GT.SAFE2 ) THEN
00350                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
00351             ELSE
00352                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
00353      $                      SAFE1
00354             END IF
00355   200    CONTINUE
00356 *
00357          KASE = 0
00358   210    CONTINUE
00359          CALL CLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE )
00360          IF( KASE.NE.0 ) THEN
00361             IF( KASE.EQ.1 ) THEN
00362 *
00363 *              Multiply by diag(W)*inv(op(A)**H).
00364 *
00365                CALL CTBSV( UPLO, TRANST, DIAG, N, KD, AB, LDAB, WORK,
00366      $                     1 )
00367                DO 220 I = 1, N
00368                   WORK( I ) = RWORK( I )*WORK( I )
00369   220          CONTINUE
00370             ELSE
00371 *
00372 *              Multiply by inv(op(A))*diag(W).
00373 *
00374                DO 230 I = 1, N
00375                   WORK( I ) = RWORK( I )*WORK( I )
00376   230          CONTINUE
00377                CALL CTBSV( UPLO, TRANSN, DIAG, N, KD, AB, LDAB, WORK,
00378      $                     1 )
00379             END IF
00380             GO TO 210
00381          END IF
00382 *
00383 *        Normalize error.
00384 *
00385          LSTRES = ZERO
00386          DO 240 I = 1, N
00387             LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) )
00388   240    CONTINUE
00389          IF( LSTRES.NE.ZERO )
00390      $      FERR( J ) = FERR( J ) / LSTRES
00391 *
00392   250 CONTINUE
00393 *
00394       RETURN
00395 *
00396 *     End of CTBRFS
00397 *
00398       END