TESTING/LIN/ctbt02.f

   1 *> \brief \b CTBT02
   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 CTBT02( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, X,
  12 *                          LDX, B, LDB, WORK, RWORK, RESID )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       CHARACTER          DIAG, TRANS, UPLO
  16 *       INTEGER            KD, LDAB, LDB, LDX, N, NRHS
  17 *       REAL               RESID
  18 *       ..
  19 *       .. Array Arguments ..
  20 *       REAL               RWORK( * )
  21 *       COMPLEX            AB( LDAB, * ), B( LDB, * ), WORK( * ),
  22 *      $                   X( LDX, * )
  23 *       ..
  24 *
  25 *
  26 *> \par Purpose:
  27 *  =============
  28 *>
  29 *> \verbatim
  30 *>
  31 *> CTBT02 computes the residual for the computed solution to a
  32 *> triangular system of linear equations  A*x = b,  A**T *x = b,  or
  33 *> A**H *x = b  when A is a triangular band matrix.  Here A**T denotes
  34 *> the transpose of A, A**H denotes the conjugate transpose of A, and
  35 *> x and b are N by NRHS matrices.  The test ratio is the maximum over
  36 *> the number of right hand sides of
  37 *>    norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
  38 *> where op(A) denotes A, A**T, or A**H, and EPS is the machine epsilon.
  39 *> \endverbatim
  40 *
  41 *  Arguments:
  42 *  ==========
  43 *
  44 *> \param[in] UPLO
  45 *> \verbatim
  46 *>          UPLO is CHARACTER*1
  47 *>          Specifies whether the matrix A is upper or lower triangular.
  48 *>          = 'U':  Upper triangular
  49 *>          = 'L':  Lower triangular
  50 *> \endverbatim
  51 *>
  52 *> \param[in] TRANS
  53 *> \verbatim
  54 *>          TRANS is CHARACTER*1
  55 *>          Specifies the operation applied to A.
  56 *>          = 'N':  A *x = b     (No transpose)
  57 *>          = 'T':  A**T *x = b  (Transpose)
  58 *>          = 'C':  A**H *x = b  (Conjugate transpose)
  59 *> \endverbatim
  60 *>
  61 *> \param[in] DIAG
  62 *> \verbatim
  63 *>          DIAG is CHARACTER*1
  64 *>          Specifies whether or not the matrix A is unit triangular.
  65 *>          = 'N':  Non-unit triangular
  66 *>          = 'U':  Unit triangular
  67 *> \endverbatim
  68 *>
  69 *> \param[in] N
  70 *> \verbatim
  71 *>          N is INTEGER
  72 *>          The order of the matrix A.  N >= 0.
  73 *> \endverbatim
  74 *>
  75 *> \param[in] KD
  76 *> \verbatim
  77 *>          KD is INTEGER
  78 *>          The number of superdiagonals or subdiagonals of the
  79 *>          triangular band matrix A.  KD >= 0.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] NRHS
  83 *> \verbatim
  84 *>          NRHS is INTEGER
  85 *>          The number of right hand sides, i.e., the number of columns
  86 *>          of the matrices X and B.  NRHS >= 0.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] AB
  90 *> \verbatim
  91 *>          AB is COMPLEX array, dimension (LDA,N)
  92 *>          The upper or lower triangular band matrix A, stored in the
  93 *>          first kd+1 rows of the array. The j-th column of A is stored
  94 *>          in the j-th column of the array AB as follows:
  95 *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
  96 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
  97 *> \endverbatim
  98 *>
  99 *> \param[in] LDAB
 100 *> \verbatim
 101 *>          LDAB is INTEGER
 102 *>          The leading dimension of the array AB.  LDAB >= max(1,KD+1).
 103 *> \endverbatim
 104 *>
 105 *> \param[in] X
 106 *> \verbatim
 107 *>          X is COMPLEX array, dimension (LDX,NRHS)
 108 *>          The computed solution vectors for the system of linear
 109 *>          equations.
 110 *> \endverbatim
 111 *>
 112 *> \param[in] LDX
 113 *> \verbatim
 114 *>          LDX is INTEGER
 115 *>          The leading dimension of the array X.  LDX >= max(1,N).
 116 *> \endverbatim
 117 *>
 118 *> \param[in] B
 119 *> \verbatim
 120 *>          B is COMPLEX array, dimension (LDB,NRHS)
 121 *>          The right hand side vectors for the system of linear
 122 *>          equations.
 123 *> \endverbatim
 124 *>
 125 *> \param[in] LDB
 126 *> \verbatim
 127 *>          LDB is INTEGER
 128 *>          The leading dimension of the array B.  LDB >= max(1,N).
 129 *> \endverbatim
 130 *>
 131 *> \param[out] WORK
 132 *> \verbatim
 133 *>          WORK is COMPLEX array, dimension (N)
 134 *> \endverbatim
 135 *>
 136 *> \param[out] RWORK
 137 *> \verbatim
 138 *>          RWORK is REAL array, dimension (N)
 139 *> \endverbatim
 140 *>
 141 *> \param[out] RESID
 142 *> \verbatim
 143 *>          RESID is REAL
 144 *>          The maximum over the number of right hand sides of
 145 *>          norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
 146 *> \endverbatim
 147 *
 148 *  Authors:
 149 *  ========
 150 *
 151 *> \author Univ. of Tennessee
 152 *> \author Univ. of California Berkeley
 153 *> \author Univ. of Colorado Denver
 154 *> \author NAG Ltd.
 155 *
 156 *> \date November 2011
 157 *
 158 *> \ingroup complex_lin
 159 *
 160 *  =====================================================================
 161       SUBROUTINE CTBT02( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, X,
 162      $                   LDX, B, LDB, WORK, RWORK, RESID )
 163 *
 164 *  -- LAPACK test routine (version 3.4.0) --
 165 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 166 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 167 *     November 2011
 168 *
 169 *     .. Scalar Arguments ..
 170       CHARACTER          DIAG, TRANS, UPLO
 171       INTEGER            KD, LDAB, LDB, LDX, N, NRHS
 172       REAL               RESID
 173 *     ..
 174 *     .. Array Arguments ..
 175       REAL               RWORK( * )
 176       COMPLEX            AB( LDAB, * ), B( LDB, * ), WORK( * ),
 177      $                   X( LDX, * )
 178 *     ..
 179 *
 180 *  =====================================================================
 181 *
 182 *     .. Parameters ..
 183       REAL               ZERO, ONE
 184       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 185 *     ..
 186 *     .. Local Scalars ..
 187       INTEGER            J
 188       REAL               ANORM, BNORM, EPS, XNORM
 189 *     ..
 190 *     .. External Functions ..
 191       LOGICAL            LSAME
 192       REAL               CLANTB, SCASUM, SLAMCH
 193       EXTERNAL           LSAME, CLANTB, SCASUM, SLAMCH
 194 *     ..
 195 *     .. External Subroutines ..
 196       EXTERNAL           CAXPY, CCOPY, CTBMV
 197 *     ..
 198 *     .. Intrinsic Functions ..
 199       INTRINSIC          CMPLX, MAX
 200 *     ..
 201 *     .. Executable Statements ..
 202 *
 203 *     Quick exit if N = 0 or NRHS = 0
 204 *
 205       IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
 206          RESID = ZERO
 207          RETURN
 208       END IF
 209 *
 210 *     Compute the 1-norm of A or A'.
 211 *
 212       IF( LSAME( TRANS, 'N' ) ) THEN
 213          ANORM = CLANTB( '1', UPLO, DIAG, N, KD, AB, LDAB, RWORK )
 214       ELSE
 215          ANORM = CLANTB( 'I', UPLO, DIAG, N, KD, AB, LDAB, RWORK )
 216       END IF
 217 *
 218 *     Exit with RESID = 1/EPS if ANORM = 0.
 219 *
 220       EPS = SLAMCH( 'Epsilon' )
 221       IF( ANORM.LE.ZERO ) THEN
 222          RESID = ONE / EPS
 223          RETURN
 224       END IF
 225 *
 226 *     Compute the maximum over the number of right hand sides of
 227 *        norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
 228 *
 229       RESID = ZERO
 230       DO 10 J = 1, NRHS
 231          CALL CCOPY( N, X( 1, J ), 1, WORK, 1 )
 232          CALL CTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK, 1 )
 233          CALL CAXPY( N, CMPLX( -ONE ), B( 1, J ), 1, WORK, 1 )
 234          BNORM = SCASUM( N, WORK, 1 )
 235          XNORM = SCASUM( N, X( 1, J ), 1 )
 236          IF( XNORM.LE.ZERO ) THEN
 237             RESID = ONE / EPS
 238          ELSE
 239             RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
 240          END IF
 241    10 CONTINUE
 242 *
 243       RETURN
 244 *
 245 *     End of CTBT02
 246 *
 247       END