TESTING/LIN/ztbt03.f

   1 *> \brief \b ZTBT03
   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 ZTBT03( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB,
  12 *                          SCALE, CNORM, TSCAL, X, LDX, B, LDB, WORK,
  13 *                          RESID )
  14 *
  15 *       .. Scalar Arguments ..
  16 *       CHARACTER          DIAG, TRANS, UPLO
  17 *       INTEGER            KD, LDAB, LDB, LDX, N, NRHS
  18 *       DOUBLE PRECISION   RESID, SCALE, TSCAL
  19 *       ..
  20 *       .. Array Arguments ..
  21 *       DOUBLE PRECISION   CNORM( * )
  22 *       COMPLEX*16         AB( LDAB, * ), B( LDB, * ), WORK( * ),
  23 *      $                   X( LDX, * )
  24 *       ..
  25 *
  26 *
  27 *> \par Purpose:
  28 *  =============
  29 *>
  30 *> \verbatim
  31 *>
  32 *> ZTBT03 computes the residual for the solution to a scaled triangular
  33 *> system of equations  A*x = s*b,  A**T *x = s*b,  or  A**H *x = s*b
  34 *> when A is a triangular band matrix.  Here A**T  denotes the transpose
  35 *> of A, A**H denotes the conjugate transpose of A, s is a scalar, and
  36 *> x and b are N by NRHS matrices.  The test ratio is the maximum over
  37 *> the number of right hand sides of
  38 *>    norm(s*b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
  39 *> where op(A) denotes A, A**T, or A**H, and EPS is the machine epsilon.
  40 *> \endverbatim
  41 *
  42 *  Arguments:
  43 *  ==========
  44 *
  45 *> \param[in] UPLO
  46 *> \verbatim
  47 *>          UPLO is CHARACTER*1
  48 *>          Specifies whether the matrix A is upper or lower triangular.
  49 *>          = 'U':  Upper triangular
  50 *>          = 'L':  Lower triangular
  51 *> \endverbatim
  52 *>
  53 *> \param[in] TRANS
  54 *> \verbatim
  55 *>          TRANS is CHARACTER*1
  56 *>          Specifies the operation applied to A.
  57 *>          = 'N':  A *x = s*b     (No transpose)
  58 *>          = 'T':  A**T *x = s*b  (Transpose)
  59 *>          = 'C':  A**H *x = s*b  (Conjugate transpose)
  60 *> \endverbatim
  61 *>
  62 *> \param[in] DIAG
  63 *> \verbatim
  64 *>          DIAG is CHARACTER*1
  65 *>          Specifies whether or not the matrix A is unit triangular.
  66 *>          = 'N':  Non-unit triangular
  67 *>          = 'U':  Unit triangular
  68 *> \endverbatim
  69 *>
  70 *> \param[in] N
  71 *> \verbatim
  72 *>          N is INTEGER
  73 *>          The order of the matrix A.  N >= 0.
  74 *> \endverbatim
  75 *>
  76 *> \param[in] KD
  77 *> \verbatim
  78 *>          KD is INTEGER
  79 *>          The number of superdiagonals or subdiagonals of the
  80 *>          triangular band matrix A.  KD >= 0.
  81 *> \endverbatim
  82 *>
  83 *> \param[in] NRHS
  84 *> \verbatim
  85 *>          NRHS is INTEGER
  86 *>          The number of right hand sides, i.e., the number of columns
  87 *>          of the matrices X and B.  NRHS >= 0.
  88 *> \endverbatim
  89 *>
  90 *> \param[in] AB
  91 *> \verbatim
  92 *>          AB is COMPLEX*16 array, dimension (LDAB,N)
  93 *>          The upper or lower triangular band matrix A, stored in the
  94 *>          first kd+1 rows of the array. The j-th column of A is stored
  95 *>          in the j-th column of the array AB as follows:
  96 *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
  97 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
  98 *> \endverbatim
  99 *>
 100 *> \param[in] LDAB
 101 *> \verbatim
 102 *>          LDAB is INTEGER
 103 *>          The leading dimension of the array AB.  LDAB >= KD+1.
 104 *> \endverbatim
 105 *>
 106 *> \param[in] SCALE
 107 *> \verbatim
 108 *>          SCALE is DOUBLE PRECISION
 109 *>          The scaling factor s used in solving the triangular system.
 110 *> \endverbatim
 111 *>
 112 *> \param[in] CNORM
 113 *> \verbatim
 114 *>          CNORM is DOUBLE PRECISION array, dimension (N)
 115 *>          The 1-norms of the columns of A, not counting the diagonal.
 116 *> \endverbatim
 117 *>
 118 *> \param[in] TSCAL
 119 *> \verbatim
 120 *>          TSCAL is DOUBLE PRECISION
 121 *>          The scaling factor used in computing the 1-norms in CNORM.
 122 *>          CNORM actually contains the column norms of TSCAL*A.
 123 *> \endverbatim
 124 *>
 125 *> \param[in] X
 126 *> \verbatim
 127 *>          X is COMPLEX*16 array, dimension (LDX,NRHS)
 128 *>          The computed solution vectors for the system of linear
 129 *>          equations.
 130 *> \endverbatim
 131 *>
 132 *> \param[in] LDX
 133 *> \verbatim
 134 *>          LDX is INTEGER
 135 *>          The leading dimension of the array X.  LDX >= max(1,N).
 136 *> \endverbatim
 137 *>
 138 *> \param[in] B
 139 *> \verbatim
 140 *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
 141 *>          The right hand side vectors for the system of linear
 142 *>          equations.
 143 *> \endverbatim
 144 *>
 145 *> \param[in] LDB
 146 *> \verbatim
 147 *>          LDB is INTEGER
 148 *>          The leading dimension of the array B.  LDB >= max(1,N).
 149 *> \endverbatim
 150 *>
 151 *> \param[out] WORK
 152 *> \verbatim
 153 *>          WORK is COMPLEX*16 array, dimension (N)
 154 *> \endverbatim
 155 *>
 156 *> \param[out] RESID
 157 *> \verbatim
 158 *>          RESID is DOUBLE PRECISION
 159 *>          The maximum over the number of right hand sides of
 160 *>          norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * 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 ZTBT03( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB,
 177      $                   SCALE, CNORM, TSCAL, X, LDX, B, LDB, WORK,
 178      $                   RESID )
 179 *
 180 *  -- LAPACK test routine (version 3.4.0) --
 181 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 182 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 183 *     November 2011
 184 *
 185 *     .. Scalar Arguments ..
 186       CHARACTER          DIAG, TRANS, UPLO
 187       INTEGER            KD, LDAB, LDB, LDX, N, NRHS
 188       DOUBLE PRECISION   RESID, SCALE, TSCAL
 189 *     ..
 190 *     .. Array Arguments ..
 191       DOUBLE PRECISION   CNORM( * )
 192       COMPLEX*16         AB( LDAB, * ), B( LDB, * ), WORK( * ),
 193      $                   X( LDX, * )
 194 *     ..
 195 *
 196 *  =====================================================================
 197 *
 198 *
 199 *     .. Parameters ..
 200       DOUBLE PRECISION   ONE, ZERO
 201       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 202 *     ..
 203 *     .. Local Scalars ..
 204       INTEGER            IX, J
 205       DOUBLE PRECISION   EPS, ERR, SMLNUM, TNORM, XNORM, XSCAL
 206 *     ..
 207 *     .. External Functions ..
 208       LOGICAL            LSAME
 209       INTEGER            IZAMAX
 210       DOUBLE PRECISION   DLAMCH
 211       EXTERNAL           LSAME, IZAMAX, DLAMCH
 212 *     ..
 213 *     .. External Subroutines ..
 214       EXTERNAL           ZAXPY, ZCOPY, ZDSCAL, ZTBMV
 215 *     ..
 216 *     .. Intrinsic Functions ..
 217       INTRINSIC          ABS, DBLE, DCMPLX, MAX
 218 *     ..
 219 *     .. Executable Statements ..
 220 *
 221 *     Quick exit if N = 0
 222 *
 223       IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
 224          RESID = ZERO
 225          RETURN
 226       END IF
 227       EPS = DLAMCH( 'Epsilon' )
 228       SMLNUM = DLAMCH( 'Safe minimum' )
 229 *
 230 *     Compute the norm of the triangular matrix A using the column
 231 *     norms already computed by ZLATBS.
 232 *
 233       TNORM = ZERO
 234       IF( LSAME( DIAG, 'N' ) ) THEN
 235          IF( LSAME( UPLO, 'U' ) ) THEN
 236             DO 10 J = 1, N
 237                TNORM = MAX( TNORM, TSCAL*ABS( AB( KD+1, J ) )+
 238      $                 CNORM( J ) )
 239    10       CONTINUE
 240          ELSE
 241             DO 20 J = 1, N
 242                TNORM = MAX( TNORM, TSCAL*ABS( AB( 1, J ) )+CNORM( J ) )
 243    20       CONTINUE
 244          END IF
 245       ELSE
 246          DO 30 J = 1, N
 247             TNORM = MAX( TNORM, TSCAL+CNORM( J ) )
 248    30    CONTINUE
 249       END IF
 250 *
 251 *     Compute the maximum over the number of right hand sides of
 252 *        norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ).
 253 *
 254       RESID = ZERO
 255       DO 40 J = 1, NRHS
 256          CALL ZCOPY( N, X( 1, J ), 1, WORK, 1 )
 257          IX = IZAMAX( N, WORK, 1 )
 258          XNORM = MAX( ONE, ABS( X( IX, J ) ) )
 259          XSCAL = ( ONE / XNORM ) / DBLE( KD+1 )
 260          CALL ZDSCAL( N, XSCAL, WORK, 1 )
 261          CALL ZTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK, 1 )
 262          CALL ZAXPY( N, DCMPLX( -SCALE*XSCAL ), B( 1, J ), 1, WORK, 1 )
 263          IX = IZAMAX( N, WORK, 1 )
 264          ERR = TSCAL*ABS( WORK( IX ) )
 265          IX = IZAMAX( N, X( 1, J ), 1 )
 266          XNORM = ABS( X( IX, J ) )
 267          IF( ERR*SMLNUM.LE.XNORM ) THEN
 268             IF( XNORM.GT.ZERO )
 269      $         ERR = ERR / XNORM
 270          ELSE
 271             IF( ERR.GT.ZERO )
 272      $         ERR = ONE / EPS
 273          END IF
 274          IF( ERR*SMLNUM.LE.TNORM ) THEN
 275             IF( TNORM.GT.ZERO )
 276      $         ERR = ERR / TNORM
 277          ELSE
 278             IF( ERR.GT.ZERO )
 279      $         ERR = ONE / EPS
 280          END IF
 281          RESID = MAX( RESID, ERR )
 282    40 CONTINUE
 283 *
 284       RETURN
 285 *
 286 *     End of ZTBT03
 287 *
 288       END