TESTING/LIN/cgtt05.f

   1 *> \brief \b CGTT05
   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 CGTT05( TRANS, N, NRHS, DL, D, DU, B, LDB, X, LDX,
  12 *                          XACT, LDXACT, FERR, BERR, RESLTS )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       CHARACTER          TRANS
  16 *       INTEGER            LDB, LDX, LDXACT, N, NRHS
  17 *       ..
  18 *       .. Array Arguments ..
  19 *       REAL               BERR( * ), FERR( * ), RESLTS( * )
  20 *       COMPLEX            B( LDB, * ), D( * ), DL( * ), DU( * ),
  21 *      $                   X( LDX, * ), XACT( LDXACT, * )
  22 *       ..
  23 *
  24 *
  25 *> \par Purpose:
  26 *  =============
  27 *>
  28 *> \verbatim
  29 *>
  30 *> CGTT05 tests the error bounds from iterative refinement for the
  31 *> computed solution to a system of equations A*X = B, where A is a
  32 *> general tridiagonal matrix of order n and op(A) = A or A**T,
  33 *> 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 and XACT.  N >= 0.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] NRHS
  65 *> \verbatim
  66 *>          NRHS is INTEGER
  67 *>          The number of columns of the matrices X and XACT.  NRHS >= 0.
  68 *> \endverbatim
  69 *>
  70 *> \param[in] DL
  71 *> \verbatim
  72 *>          DL is COMPLEX array, dimension (N-1)
  73 *>          The (n-1) sub-diagonal elements of A.
  74 *> \endverbatim
  75 *>
  76 *> \param[in] D
  77 *> \verbatim
  78 *>          D is COMPLEX array, dimension (N)
  79 *>          The diagonal elements of A.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] DU
  83 *> \verbatim
  84 *>          DU is COMPLEX array, dimension (N-1)
  85 *>          The (n-1) super-diagonal elements of A.
  86 *> \endverbatim
  87 *>
  88 *> \param[in] B
  89 *> \verbatim
  90 *>          B is COMPLEX array, dimension (LDB,NRHS)
  91 *>          The right hand side vectors for the system of linear
  92 *>          equations.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] LDB
  96 *> \verbatim
  97 *>          LDB is INTEGER
  98 *>          The leading dimension of the array B.  LDB >= max(1,N).
  99 *> \endverbatim
 100 *>
 101 *> \param[in] X
 102 *> \verbatim
 103 *>          X is COMPLEX array, dimension (LDX,NRHS)
 104 *>          The computed solution vectors.  Each vector is stored as a
 105 *>          column of the matrix X.
 106 *> \endverbatim
 107 *>
 108 *> \param[in] LDX
 109 *> \verbatim
 110 *>          LDX is INTEGER
 111 *>          The leading dimension of the array X.  LDX >= max(1,N).
 112 *> \endverbatim
 113 *>
 114 *> \param[in] XACT
 115 *> \verbatim
 116 *>          XACT is COMPLEX array, dimension (LDX,NRHS)
 117 *>          The exact solution vectors.  Each vector is stored as a
 118 *>          column of the matrix XACT.
 119 *> \endverbatim
 120 *>
 121 *> \param[in] LDXACT
 122 *> \verbatim
 123 *>          LDXACT is INTEGER
 124 *>          The leading dimension of the array XACT.  LDXACT >= max(1,N).
 125 *> \endverbatim
 126 *>
 127 *> \param[in] FERR
 128 *> \verbatim
 129 *>          FERR is REAL array, dimension (NRHS)
 130 *>          The estimated forward error bounds for each solution vector
 131 *>          X.  If XTRUE is the true solution, FERR bounds the magnitude
 132 *>          of the largest entry in (X - XTRUE) divided by the magnitude
 133 *>          of the largest entry in X.
 134 *> \endverbatim
 135 *>
 136 *> \param[in] BERR
 137 *> \verbatim
 138 *>          BERR is REAL 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 REAL 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 / ( NZ*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 complex_lin
 163 *
 164 *  =====================================================================
 165       SUBROUTINE CGTT05( TRANS, N, NRHS, DL, D, DU, B, LDB, X, LDX,
 166      $                   XACT, LDXACT, FERR, 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       INTEGER            LDB, LDX, LDXACT, N, NRHS
 176 *     ..
 177 *     .. Array Arguments ..
 178       REAL               BERR( * ), FERR( * ), RESLTS( * )
 179       COMPLEX            B( LDB, * ), D( * ), DL( * ), DU( * ),
 180      $                   X( LDX, * ), XACT( LDXACT, * )
 181 *     ..
 182 *
 183 *  =====================================================================
 184 *
 185 *     .. Parameters ..
 186       REAL               ZERO, ONE
 187       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 188 *     ..
 189 *     .. Local Scalars ..
 190       LOGICAL            NOTRAN
 191       INTEGER            I, IMAX, J, K, NZ
 192       REAL               AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
 193       COMPLEX            ZDUM
 194 *     ..
 195 *     .. External Functions ..
 196       LOGICAL            LSAME
 197       INTEGER            ICAMAX
 198       REAL               SLAMCH
 199       EXTERNAL           LSAME, ICAMAX, SLAMCH
 200 *     ..
 201 *     .. Intrinsic Functions ..
 202       INTRINSIC          ABS, AIMAG, MAX, MIN, REAL
 203 *     ..
 204 *     .. Statement Functions ..
 205       REAL               CABS1
 206 *     ..
 207 *     .. Statement Function definitions ..
 208       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
 209 *     ..
 210 *     .. Executable Statements ..
 211 *
 212 *     Quick exit if N = 0 or NRHS = 0.
 213 *
 214       IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
 215          RESLTS( 1 ) = ZERO
 216          RESLTS( 2 ) = ZERO
 217          RETURN
 218       END IF
 219 *
 220       EPS = SLAMCH( 'Epsilon' )
 221       UNFL = SLAMCH( 'Safe minimum' )
 222       OVFL = ONE / UNFL
 223       NOTRAN = LSAME( TRANS, 'N' )
 224       NZ = 4
 225 *
 226 *     Test 1:  Compute the maximum of
 227 *        norm(X - XACT) / ( norm(X) * FERR )
 228 *     over all the vectors X and XACT using the infinity-norm.
 229 *
 230       ERRBND = ZERO
 231       DO 30 J = 1, NRHS
 232          IMAX = ICAMAX( N, X( 1, J ), 1 )
 233          XNORM = MAX( CABS1( X( IMAX, J ) ), UNFL )
 234          DIFF = ZERO
 235          DO 10 I = 1, N
 236             DIFF = MAX( DIFF, CABS1( X( I, J )-XACT( I, J ) ) )
 237    10    CONTINUE
 238 *
 239          IF( XNORM.GT.ONE ) THEN
 240             GO TO 20
 241          ELSE IF( DIFF.LE.OVFL*XNORM ) THEN
 242             GO TO 20
 243          ELSE
 244             ERRBND = ONE / EPS
 245             GO TO 30
 246          END IF
 247 *
 248    20    CONTINUE
 249          IF( DIFF / XNORM.LE.FERR( J ) ) THEN
 250             ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) )
 251          ELSE
 252             ERRBND = ONE / EPS
 253          END IF
 254    30 CONTINUE
 255       RESLTS( 1 ) = ERRBND
 256 *
 257 *     Test 2:  Compute the maximum of BERR / ( NZ*EPS + (*) ), where
 258 *     (*) = NZ*UNFL / (min_i (abs(op(A))*abs(X) +abs(b))_i )
 259 *
 260       DO 60 K = 1, NRHS
 261          IF( NOTRAN ) THEN
 262             IF( N.EQ.1 ) THEN
 263                AXBI = CABS1( B( 1, K ) ) +
 264      $                CABS1( D( 1 ) )*CABS1( X( 1, K ) )
 265             ELSE
 266                AXBI = CABS1( B( 1, K ) ) +
 267      $                CABS1( D( 1 ) )*CABS1( X( 1, K ) ) +
 268      $                CABS1( DU( 1 ) )*CABS1( X( 2, K ) )
 269                DO 40 I = 2, N - 1
 270                   TMP = CABS1( B( I, K ) ) +
 271      $                  CABS1( DL( I-1 ) )*CABS1( X( I-1, K ) ) +
 272      $                  CABS1( D( I ) )*CABS1( X( I, K ) ) +
 273      $                  CABS1( DU( I ) )*CABS1( X( I+1, K ) )
 274                   AXBI = MIN( AXBI, TMP )
 275    40          CONTINUE
 276                TMP = CABS1( B( N, K ) ) + CABS1( DL( N-1 ) )*
 277      $               CABS1( X( N-1, K ) ) + CABS1( D( N ) )*
 278      $               CABS1( X( N, K ) )
 279                AXBI = MIN( AXBI, TMP )
 280             END IF
 281          ELSE
 282             IF( N.EQ.1 ) THEN
 283                AXBI = CABS1( B( 1, K ) ) +
 284      $                CABS1( D( 1 ) )*CABS1( X( 1, K ) )
 285             ELSE
 286                AXBI = CABS1( B( 1, K ) ) +
 287      $                CABS1( D( 1 ) )*CABS1( X( 1, K ) ) +
 288      $                CABS1( DL( 1 ) )*CABS1( X( 2, K ) )
 289                DO 50 I = 2, N - 1
 290                   TMP = CABS1( B( I, K ) ) +
 291      $                  CABS1( DU( I-1 ) )*CABS1( X( I-1, K ) ) +
 292      $                  CABS1( D( I ) )*CABS1( X( I, K ) ) +
 293      $                  CABS1( DL( I ) )*CABS1( X( I+1, K ) )
 294                   AXBI = MIN( AXBI, TMP )
 295    50          CONTINUE
 296                TMP = CABS1( B( N, K ) ) + CABS1( DU( N-1 ) )*
 297      $               CABS1( X( N-1, K ) ) + CABS1( D( N ) )*
 298      $               CABS1( X( N, K ) )
 299                AXBI = MIN( AXBI, TMP )
 300             END IF
 301          END IF
 302          TMP = BERR( K ) / ( NZ*EPS+NZ*UNFL / MAX( AXBI, NZ*UNFL ) )
 303          IF( K.EQ.1 ) THEN
 304             RESLTS( 2 ) = TMP
 305          ELSE
 306             RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP )
 307          END IF
 308    60 CONTINUE
 309 *
 310       RETURN
 311 *
 312 *     End of CGTT05
 313 *
 314       END