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