TESTING/LIN/zgbt05.f

   1 *> \brief \b ZGBT05
   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 ZGBT05( 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 *       DOUBLE PRECISION   BERR( * ), FERR( * ), RESLTS( * )
  20 *       COMPLEX*16         AB( LDAB, * ), B( LDB, * ), X( LDX, * ),
  21 *      $                   XACT( LDXACT, * )
  22 *       ..
  23 *
  24 *
  25 *> \par Purpose:
  26 *  =============
  27 *>
  28 *> \verbatim
  29 *>
  30 *> ZGBT05 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 complex16_lin
 174 *
 175 *  =====================================================================
 176       SUBROUTINE ZGBT05( 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       DOUBLE PRECISION   BERR( * ), FERR( * ), RESLTS( * )
 190       COMPLEX*16         AB( LDAB, * ), B( LDB, * ), X( LDX, * ),
 191      $                   XACT( LDXACT, * )
 192 *     ..
 193 *
 194 *  =====================================================================
 195 *
 196 *     .. Parameters ..
 197       DOUBLE PRECISION   ZERO, ONE
 198       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 199 *     ..
 200 *     .. Local Scalars ..
 201       LOGICAL            NOTRAN
 202       INTEGER            I, IMAX, J, K, NZ
 203       DOUBLE PRECISION   AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
 204       COMPLEX*16         ZDUM
 205 *     ..
 206 *     .. External Functions ..
 207       LOGICAL            LSAME
 208       INTEGER            IZAMAX
 209       DOUBLE PRECISION   DLAMCH
 210       EXTERNAL           LSAME, IZAMAX, DLAMCH
 211 *     ..
 212 *     .. Intrinsic Functions ..
 213       INTRINSIC          ABS, DBLE, DIMAG, MAX, MIN
 214 *     ..
 215 *     .. Statement Functions ..
 216       DOUBLE PRECISION   CABS1
 217 *     ..
 218 *     .. Statement Function definitions ..
 219       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
 220 *     ..
 221 *     .. Executable Statements ..
 222 *
 223 *     Quick exit if N = 0 or NRHS = 0.
 224 *
 225       IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
 226          RESLTS( 1 ) = ZERO
 227          RESLTS( 2 ) = ZERO
 228          RETURN
 229       END IF
 230 *
 231       EPS = DLAMCH( 'Epsilon' )
 232       UNFL = DLAMCH( 'Safe minimum' )
 233       OVFL = ONE / UNFL
 234       NOTRAN = LSAME( TRANS, 'N' )
 235       NZ = MIN( KL+KU+2, N+1 )
 236 *
 237 *     Test 1:  Compute the maximum of
 238 *        norm(X - XACT) / ( norm(X) * FERR )
 239 *     over all the vectors X and XACT using the infinity-norm.
 240 *
 241       ERRBND = ZERO
 242       DO 30 J = 1, NRHS
 243          IMAX = IZAMAX( N, X( 1, J ), 1 )
 244          XNORM = MAX( CABS1( X( IMAX, J ) ), UNFL )
 245          DIFF = ZERO
 246          DO 10 I = 1, N
 247             DIFF = MAX( DIFF, CABS1( X( I, J )-XACT( I, J ) ) )
 248    10    CONTINUE
 249 *
 250          IF( XNORM.GT.ONE ) THEN
 251             GO TO 20
 252          ELSE IF( DIFF.LE.OVFL*XNORM ) THEN
 253             GO TO 20
 254          ELSE
 255             ERRBND = ONE / EPS
 256             GO TO 30
 257          END IF
 258 *
 259    20    CONTINUE
 260          IF( DIFF / XNORM.LE.FERR( J ) ) THEN
 261             ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) )
 262          ELSE
 263             ERRBND = ONE / EPS
 264          END IF
 265    30 CONTINUE
 266       RESLTS( 1 ) = ERRBND
 267 *
 268 *     Test 2:  Compute the maximum of BERR / ( NZ*EPS + (*) ), where
 269 *     (*) = NZ*UNFL / (min_i (abs(op(A))*abs(X) +abs(b))_i )
 270 *
 271       DO 70 K = 1, NRHS
 272          DO 60 I = 1, N
 273             TMP = CABS1( B( I, K ) )
 274             IF( NOTRAN ) THEN
 275                DO 40 J = MAX( I-KL, 1 ), MIN( I+KU, N )
 276                   TMP = TMP + CABS1( AB( KU+1+I-J, J ) )*
 277      $                  CABS1( X( J, K ) )
 278    40          CONTINUE
 279             ELSE
 280                DO 50 J = MAX( I-KU, 1 ), MIN( I+KL, N )
 281                   TMP = TMP + CABS1( AB( KU+1+J-I, I ) )*
 282      $                  CABS1( X( J, K ) )
 283    50          CONTINUE
 284             END IF
 285             IF( I.EQ.1 ) THEN
 286                AXBI = TMP
 287             ELSE
 288                AXBI = MIN( AXBI, TMP )
 289             END IF
 290    60    CONTINUE
 291          TMP = BERR( K ) / ( NZ*EPS+NZ*UNFL / MAX( AXBI, NZ*UNFL ) )
 292          IF( K.EQ.1 ) THEN
 293             RESLTS( 2 ) = TMP
 294          ELSE
 295             RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP )
 296          END IF
 297    70 CONTINUE
 298 *
 299       RETURN
 300 *
 301 *     End of ZGBT05
 302 *
 303       END