TESTING/LIN/dget07.f

   1 *> \brief \b DGET07
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE DGET07( TRANS, N, NRHS, A, LDA, B, LDB, X, LDX, XACT,
  12 *                          LDXACT, FERR, CHKFERR, BERR, RESLTS )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       CHARACTER          TRANS
  16 *       LOGICAL            CHKFERR
  17 *       INTEGER            LDA, LDB, LDX, LDXACT, N, NRHS
  18 *       ..
  19 *       .. Array Arguments ..
  20 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ),
  21 *      $                   RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
  22 *       ..
  23 *
  24 *
  25 *> \par Purpose:
  26 *  =============
  27 *>
  28 *> \verbatim
  29 *>
  30 *> DGET07 tests the error bounds from iterative refinement for the
  31 *> computed solution to a system of equations op(A)*X = B, where A is a
  32 *> general n by n matrix and op(A) = A or A**T, depending on TRANS.
  33 *>
  34 *> RESLTS(1) = test of the error bound
  35 *>           = norm(X - XACT) / ( norm(X) * FERR )
  36 *>
  37 *> A large value is returned if this ratio is not less than one.
  38 *>
  39 *> RESLTS(2) = residual from the iterative refinement routine
  40 *>           = the maximum of BERR / ( (n+1)*EPS + (*) ), where
  41 *>             (*) = (n+1)*UNFL / (min_i (abs(op(A))*abs(X) +abs(b))_i )
  42 *> \endverbatim
  43 *
  44 *  Arguments:
  45 *  ==========
  46 *
  47 *> \param[in] TRANS
  48 *> \verbatim
  49 *>          TRANS is CHARACTER*1
  50 *>          Specifies the form of the system of equations.
  51 *>          = 'N':  A * X = B     (No transpose)
  52 *>          = 'T':  A**T * X = B  (Transpose)
  53 *>          = 'C':  A**H * X = B  (Conjugate transpose = Transpose)
  54 *> \endverbatim
  55 *>
  56 *> \param[in] N
  57 *> \verbatim
  58 *>          N is INTEGER
  59 *>          The number of rows of the matrices X and XACT.  N >= 0.
  60 *> \endverbatim
  61 *>
  62 *> \param[in] NRHS
  63 *> \verbatim
  64 *>          NRHS is INTEGER
  65 *>          The number of columns of the matrices X and XACT.  NRHS >= 0.
  66 *> \endverbatim
  67 *>
  68 *> \param[in] A
  69 *> \verbatim
  70 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
  71 *>          The original n by n matrix A.
  72 *> \endverbatim
  73 *>
  74 *> \param[in] LDA
  75 *> \verbatim
  76 *>          LDA is INTEGER
  77 *>          The leading dimension of the array A.  LDA >= max(1,N).
  78 *> \endverbatim
  79 *>
  80 *> \param[in] B
  81 *> \verbatim
  82 *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
  83 *>          The right hand side vectors for the system of linear
  84 *>          equations.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] LDB
  88 *> \verbatim
  89 *>          LDB is INTEGER
  90 *>          The leading dimension of the array B.  LDB >= max(1,N).
  91 *> \endverbatim
  92 *>
  93 *> \param[in] X
  94 *> \verbatim
  95 *>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
  96 *>          The computed solution vectors.  Each vector is stored as a
  97 *>          column of the matrix X.
  98 *> \endverbatim
  99 *>
 100 *> \param[in] LDX
 101 *> \verbatim
 102 *>          LDX is INTEGER
 103 *>          The leading dimension of the array X.  LDX >= max(1,N).
 104 *> \endverbatim
 105 *>
 106 *> \param[in] XACT
 107 *> \verbatim
 108 *>          XACT is DOUBLE PRECISION array, dimension (LDX,NRHS)
 109 *>          The exact solution vectors.  Each vector is stored as a
 110 *>          column of the matrix XACT.
 111 *> \endverbatim
 112 *>
 113 *> \param[in] LDXACT
 114 *> \verbatim
 115 *>          LDXACT is INTEGER
 116 *>          The leading dimension of the array XACT.  LDXACT >= max(1,N).
 117 *> \endverbatim
 118 *>
 119 *> \param[in] FERR
 120 *> \verbatim
 121 *>          FERR is DOUBLE PRECISION array, dimension (NRHS)
 122 *>          The estimated forward error bounds for each solution vector
 123 *>          X.  If XTRUE is the true solution, FERR bounds the magnitude
 124 *>          of the largest entry in (X - XTRUE) divided by the magnitude
 125 *>          of the largest entry in X.
 126 *> \endverbatim
 127 *>
 128 *> \param[in] CHKFERR
 129 *> \verbatim
 130 *>          CHKFERR is LOGICAL
 131 *>          Set to .TRUE. to check FERR, .FALSE. not to check FERR.
 132 *>          When the test system is ill-conditioned, the "true"
 133 *>          solution in XACT may be incorrect.
 134 *> \endverbatim
 135 *>
 136 *> \param[in] BERR
 137 *> \verbatim
 138 *>          BERR is DOUBLE PRECISION array, dimension (NRHS)
 139 *>          The componentwise relative backward error of each solution
 140 *>          vector (i.e., the smallest relative change in any entry of A
 141 *>          or B that makes X an exact solution).
 142 *> \endverbatim
 143 *>
 144 *> \param[out] RESLTS
 145 *> \verbatim
 146 *>          RESLTS is DOUBLE PRECISION array, dimension (2)
 147 *>          The maximum over the NRHS solution vectors of the ratios:
 148 *>          RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
 149 *>          RESLTS(2) = BERR / ( (n+1)*EPS + (*) )
 150 *> \endverbatim
 151 *
 152 *  Authors:
 153 *  ========
 154 *
 155 *> \author Univ. of Tennessee
 156 *> \author Univ. of California Berkeley
 157 *> \author Univ. of Colorado Denver
 158 *> \author NAG Ltd.
 159 *
 160 *> \date November 2011
 161 *
 162 *> \ingroup double_lin
 163 *
 164 *  =====================================================================
 165       SUBROUTINE DGET07( TRANS, N, NRHS, A, LDA, B, LDB, X, LDX, XACT,
 166      $                   LDXACT, FERR, CHKFERR, BERR, RESLTS )
 167 *
 168 *  -- LAPACK test routine (version 3.4.0) --
 169 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 170 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 171 *     November 2011
 172 *
 173 *     .. Scalar Arguments ..
 174       CHARACTER          TRANS
 175       LOGICAL            CHKFERR
 176       INTEGER            LDA, LDB, LDX, LDXACT, N, NRHS
 177 *     ..
 178 *     .. Array Arguments ..
 179       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ),
 180      $                   RESLTS( * ), X( LDX, * ), XACT( LDXACT, * )
 181 *     ..
 182 *
 183 *  =====================================================================
 184 *
 185 *     .. Parameters ..
 186       DOUBLE PRECISION   ZERO, ONE
 187       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 188 *     ..
 189 *     .. Local Scalars ..
 190       LOGICAL            NOTRAN
 191       INTEGER            I, IMAX, J, K
 192       DOUBLE PRECISION   AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
 193 *     ..
 194 *     .. External Functions ..
 195       LOGICAL            LSAME
 196       INTEGER            IDAMAX
 197       DOUBLE PRECISION   DLAMCH
 198       EXTERNAL           LSAME, IDAMAX, DLAMCH
 199 *     ..
 200 *     .. Intrinsic Functions ..
 201       INTRINSIC          ABS, MAX, MIN
 202 *     ..
 203 *     .. Executable Statements ..
 204 *
 205 *     Quick exit if N = 0 or NRHS = 0.
 206 *
 207       IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
 208          RESLTS( 1 ) = ZERO
 209          RESLTS( 2 ) = ZERO
 210          RETURN
 211       END IF
 212 *
 213       EPS = DLAMCH( 'Epsilon' )
 214       UNFL = DLAMCH( 'Safe minimum' )
 215       OVFL = ONE / UNFL
 216       NOTRAN = LSAME( TRANS, 'N' )
 217 *
 218 *     Test 1:  Compute the maximum of
 219 *        norm(X - XACT) / ( norm(X) * FERR )
 220 *     over all the vectors X and XACT using the infinity-norm.
 221 *
 222       ERRBND = ZERO
 223       IF( CHKFERR ) THEN
 224          DO 30 J = 1, NRHS
 225             IMAX = IDAMAX( N, X( 1, J ), 1 )
 226             XNORM = MAX( ABS( X( IMAX, J ) ), UNFL )
 227             DIFF = ZERO
 228             DO 10 I = 1, N
 229                DIFF = MAX( DIFF, ABS( X( I, J )-XACT( I, J ) ) )
 230  10         CONTINUE
 231 *
 232             IF( XNORM.GT.ONE ) THEN
 233                GO TO 20
 234             ELSE IF( DIFF.LE.OVFL*XNORM ) THEN
 235                GO TO 20
 236             ELSE
 237                ERRBND = ONE / EPS
 238                GO TO 30
 239             END IF
 240 *
 241  20         CONTINUE
 242             IF( DIFF / XNORM.LE.FERR( J ) ) THEN
 243                ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) )
 244             ELSE
 245                ERRBND = ONE / EPS
 246             END IF
 247  30      CONTINUE
 248       END IF
 249       RESLTS( 1 ) = ERRBND
 250 *
 251 *     Test 2:  Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where
 252 *     (*) = (n+1)*UNFL / (min_i (abs(op(A))*abs(X) +abs(b))_i )
 253 *
 254       DO 70 K = 1, NRHS
 255          DO 60 I = 1, N
 256             TMP = ABS( B( I, K ) )
 257             IF( NOTRAN ) THEN
 258                DO 40 J = 1, N
 259                   TMP = TMP + ABS( A( I, J ) )*ABS( X( J, K ) )
 260    40          CONTINUE
 261             ELSE
 262                DO 50 J = 1, N
 263                   TMP = TMP + ABS( A( J, I ) )*ABS( X( J, K ) )
 264    50          CONTINUE
 265             END IF
 266             IF( I.EQ.1 ) THEN
 267                AXBI = TMP
 268             ELSE
 269                AXBI = MIN( AXBI, TMP )
 270             END IF
 271    60    CONTINUE
 272          TMP = BERR( K ) / ( ( N+1 )*EPS+( N+1 )*UNFL /
 273      $         MAX( AXBI, ( N+1 )*UNFL ) )
 274          IF( K.EQ.1 ) THEN
 275             RESLTS( 2 ) = TMP
 276          ELSE
 277             RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP )
 278          END IF
 279    70 CONTINUE
 280 *
 281       RETURN
 282 *
 283 *     End of DGET07
 284 *
 285       END