LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dtrt05.f
Go to the documentation of this file.
00001 *> \brief \b DTRT05
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE DTRT05( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
00012 *                          LDX, XACT, LDXACT, FERR, BERR, RESLTS )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          DIAG, TRANS, UPLO
00016 *       INTEGER            LDA, LDB, LDX, LDXACT, N, NRHS
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ),
00020 *      $                   RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> DTRT05 tests the error bounds from iterative refinement for the
00030 *> computed solution to a system of equations A*X = B, where A is a
00031 *> triangular n by n matrix.
00032 *>
00033 *> RESLTS(1) = test of the error bound
00034 *>           = norm(X - XACT) / ( norm(X) * FERR )
00035 *>
00036 *> A large value is returned if this ratio is not less than one.
00037 *>
00038 *> RESLTS(2) = residual from the iterative refinement routine
00039 *>           = the maximum of BERR / ( (n+1)*EPS + (*) ), where
00040 *>             (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
00041 *> \endverbatim
00042 *
00043 *  Arguments:
00044 *  ==========
00045 *
00046 *> \param[in] UPLO
00047 *> \verbatim
00048 *>          UPLO is CHARACTER*1
00049 *>          Specifies whether the matrix A is upper or lower triangular.
00050 *>          = 'U':  Upper triangular
00051 *>          = 'L':  Lower triangular
00052 *> \endverbatim
00053 *>
00054 *> \param[in] TRANS
00055 *> \verbatim
00056 *>          TRANS is CHARACTER*1
00057 *>          Specifies the form of the system of equations.
00058 *>          = 'N':  A * X = B  (No transpose)
00059 *>          = 'T':  A'* X = B  (Transpose)
00060 *>          = 'C':  A'* X = B  (Conjugate transpose = Transpose)
00061 *> \endverbatim
00062 *>
00063 *> \param[in] DIAG
00064 *> \verbatim
00065 *>          DIAG is CHARACTER*1
00066 *>          Specifies whether or not the matrix A is unit triangular.
00067 *>          = 'N':  Non-unit triangular
00068 *>          = 'U':  Unit triangular
00069 *> \endverbatim
00070 *>
00071 *> \param[in] N
00072 *> \verbatim
00073 *>          N is INTEGER
00074 *>          The number of rows of the matrices X, B, and XACT, and the
00075 *>          order of the matrix A.  N >= 0.
00076 *> \endverbatim
00077 *>
00078 *> \param[in] NRHS
00079 *> \verbatim
00080 *>          NRHS is INTEGER
00081 *>          The number of columns of the matrices X, B, and XACT.
00082 *>          NRHS >= 0.
00083 *> \endverbatim
00084 *>
00085 *> \param[in] A
00086 *> \verbatim
00087 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
00088 *>          The triangular matrix A.  If UPLO = 'U', the leading n by n
00089 *>          upper triangular part of the array A contains the upper
00090 *>          triangular matrix, and the strictly lower triangular part of
00091 *>          A is not referenced.  If UPLO = 'L', the leading n by n lower
00092 *>          triangular part of the array A contains the lower triangular
00093 *>          matrix, and the strictly upper triangular part of A is not
00094 *>          referenced.  If DIAG = 'U', the diagonal elements of A are
00095 *>          also not referenced and are assumed to be 1.
00096 *> \endverbatim
00097 *>
00098 *> \param[in] LDA
00099 *> \verbatim
00100 *>          LDA is INTEGER
00101 *>          The leading dimension of the array A.  LDA >= max(1,N).
00102 *> \endverbatim
00103 *>
00104 *> \param[in] B
00105 *> \verbatim
00106 *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
00107 *>          The right hand side vectors for the system of linear
00108 *>          equations.
00109 *> \endverbatim
00110 *>
00111 *> \param[in] LDB
00112 *> \verbatim
00113 *>          LDB is INTEGER
00114 *>          The leading dimension of the array B.  LDB >= max(1,N).
00115 *> \endverbatim
00116 *>
00117 *> \param[in] X
00118 *> \verbatim
00119 *>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
00120 *>          The computed solution vectors.  Each vector is stored as a
00121 *>          column of the matrix X.
00122 *> \endverbatim
00123 *>
00124 *> \param[in] LDX
00125 *> \verbatim
00126 *>          LDX is INTEGER
00127 *>          The leading dimension of the array X.  LDX >= max(1,N).
00128 *> \endverbatim
00129 *>
00130 *> \param[in] XACT
00131 *> \verbatim
00132 *>          XACT is DOUBLE PRECISION array, dimension (LDX,NRHS)
00133 *>          The exact solution vectors.  Each vector is stored as a
00134 *>          column of the matrix XACT.
00135 *> \endverbatim
00136 *>
00137 *> \param[in] LDXACT
00138 *> \verbatim
00139 *>          LDXACT is INTEGER
00140 *>          The leading dimension of the array XACT.  LDXACT >= max(1,N).
00141 *> \endverbatim
00142 *>
00143 *> \param[in] FERR
00144 *> \verbatim
00145 *>          FERR is DOUBLE PRECISION array, dimension (NRHS)
00146 *>          The estimated forward error bounds for each solution vector
00147 *>          X.  If XTRUE is the true solution, FERR bounds the magnitude
00148 *>          of the largest entry in (X - XTRUE) divided by the magnitude
00149 *>          of the largest entry in X.
00150 *> \endverbatim
00151 *>
00152 *> \param[in] BERR
00153 *> \verbatim
00154 *>          BERR is DOUBLE PRECISION array, dimension (NRHS)
00155 *>          The componentwise relative backward error of each solution
00156 *>          vector (i.e., the smallest relative change in any entry of A
00157 *>          or B that makes X an exact solution).
00158 *> \endverbatim
00159 *>
00160 *> \param[out] RESLTS
00161 *> \verbatim
00162 *>          RESLTS is DOUBLE PRECISION array, dimension (2)
00163 *>          The maximum over the NRHS solution vectors of the ratios:
00164 *>          RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
00165 *>          RESLTS(2) = BERR / ( (n+1)*EPS + (*) )
00166 *> \endverbatim
00167 *
00168 *  Authors:
00169 *  ========
00170 *
00171 *> \author Univ. of Tennessee 
00172 *> \author Univ. of California Berkeley 
00173 *> \author Univ. of Colorado Denver 
00174 *> \author NAG Ltd. 
00175 *
00176 *> \date November 2011
00177 *
00178 *> \ingroup double_lin
00179 *
00180 *  =====================================================================
00181       SUBROUTINE DTRT05( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
00182      $                   LDX, XACT, LDXACT, FERR, BERR, RESLTS )
00183 *
00184 *  -- LAPACK test routine (version 3.4.0) --
00185 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00186 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00187 *     November 2011
00188 *
00189 *     .. Scalar Arguments ..
00190       CHARACTER          DIAG, TRANS, UPLO
00191       INTEGER            LDA, LDB, LDX, LDXACT, N, NRHS
00192 *     ..
00193 *     .. Array Arguments ..
00194       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ),
00195      $                   RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
00196 *     ..
00197 *
00198 *  =====================================================================
00199 *
00200 *     .. Parameters ..
00201       DOUBLE PRECISION   ZERO, ONE
00202       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00203 *     ..
00204 *     .. Local Scalars ..
00205       LOGICAL            NOTRAN, UNIT, UPPER
00206       INTEGER            I, IFU, IMAX, J, K
00207       DOUBLE PRECISION   AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
00208 *     ..
00209 *     .. External Functions ..
00210       LOGICAL            LSAME
00211       INTEGER            IDAMAX
00212       DOUBLE PRECISION   DLAMCH
00213       EXTERNAL           LSAME, IDAMAX, DLAMCH
00214 *     ..
00215 *     .. Intrinsic Functions ..
00216       INTRINSIC          ABS, MAX, MIN
00217 *     ..
00218 *     .. Executable Statements ..
00219 *
00220 *     Quick exit if N = 0 or NRHS = 0.
00221 *
00222       IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
00223          RESLTS( 1 ) = ZERO
00224          RESLTS( 2 ) = ZERO
00225          RETURN
00226       END IF
00227 *
00228       EPS = DLAMCH( 'Epsilon' )
00229       UNFL = DLAMCH( 'Safe minimum' )
00230       OVFL = ONE / UNFL
00231       UPPER = LSAME( UPLO, 'U' )
00232       NOTRAN = LSAME( TRANS, 'N' )
00233       UNIT = LSAME( DIAG, 'U' )
00234 *
00235 *     Test 1:  Compute the maximum of
00236 *        norm(X - XACT) / ( norm(X) * FERR )
00237 *     over all the vectors X and XACT using the infinity-norm.
00238 *
00239       ERRBND = ZERO
00240       DO 30 J = 1, NRHS
00241          IMAX = IDAMAX( N, X( 1, J ), 1 )
00242          XNORM = MAX( ABS( X( IMAX, J ) ), UNFL )
00243          DIFF = ZERO
00244          DO 10 I = 1, N
00245             DIFF = MAX( DIFF, ABS( X( I, J )-XACT( I, J ) ) )
00246    10    CONTINUE
00247 *
00248          IF( XNORM.GT.ONE ) THEN
00249             GO TO 20
00250          ELSE IF( DIFF.LE.OVFL*XNORM ) THEN
00251             GO TO 20
00252          ELSE
00253             ERRBND = ONE / EPS
00254             GO TO 30
00255          END IF
00256 *
00257    20    CONTINUE
00258          IF( DIFF / XNORM.LE.FERR( J ) ) THEN
00259             ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) )
00260          ELSE
00261             ERRBND = ONE / EPS
00262          END IF
00263    30 CONTINUE
00264       RESLTS( 1 ) = ERRBND
00265 *
00266 *     Test 2:  Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
00267 *     (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
00268 *
00269       IFU = 0
00270       IF( UNIT )
00271      $   IFU = 1
00272       DO 90 K = 1, NRHS
00273          DO 80 I = 1, N
00274             TMP = ABS( B( I, K ) )
00275             IF( UPPER ) THEN
00276                IF( .NOT.NOTRAN ) THEN
00277                   DO 40 J = 1, I - IFU
00278                      TMP = TMP + ABS( A( J, I ) )*ABS( X( J, K ) )
00279    40             CONTINUE
00280                   IF( UNIT )
00281      $               TMP = TMP + ABS( X( I, K ) )
00282                ELSE
00283                   IF( UNIT )
00284      $               TMP = TMP + ABS( X( I, K ) )
00285                   DO 50 J = I + IFU, N
00286                      TMP = TMP + ABS( A( I, J ) )*ABS( X( J, K ) )
00287    50             CONTINUE
00288                END IF
00289             ELSE
00290                IF( NOTRAN ) THEN
00291                   DO 60 J = 1, I - IFU
00292                      TMP = TMP + ABS( A( I, J ) )*ABS( X( J, K ) )
00293    60             CONTINUE
00294                   IF( UNIT )
00295      $               TMP = TMP + ABS( X( I, K ) )
00296                ELSE
00297                   IF( UNIT )
00298      $               TMP = TMP + ABS( X( I, K ) )
00299                   DO 70 J = I + IFU, N
00300                      TMP = TMP + ABS( A( J, I ) )*ABS( X( J, K ) )
00301    70             CONTINUE
00302                END IF
00303             END IF
00304             IF( I.EQ.1 ) THEN
00305                AXBI = TMP
00306             ELSE
00307                AXBI = MIN( AXBI, TMP )
00308             END IF
00309    80    CONTINUE
00310          TMP = BERR( K ) / ( ( N+1 )*EPS+( N+1 )*UNFL /
00311      $         MAX( AXBI, ( N+1 )*UNFL ) )
00312          IF( K.EQ.1 ) THEN
00313             RESLTS( 2 ) = TMP
00314          ELSE
00315             RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP )
00316          END IF
00317    90 CONTINUE
00318 *
00319       RETURN
00320 *
00321 *     End of DTRT05
00322 *
00323       END
 All Files Functions