TESTING/LIN/sgbt05.f

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