TESTING/LIN/ctrt02.f

   1 *> \brief \b CTRT02
   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 CTRT02( UPLO, TRANS, DIAG, N, NRHS, A, LDA, X, LDX, B,
  12 *                          LDB, WORK, RWORK, RESID )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       CHARACTER          DIAG, TRANS, UPLO
  16 *       INTEGER            LDA, LDB, LDX, N, NRHS
  17 *       REAL               RESID
  18 *       ..
  19 *       .. Array Arguments ..
  20 *       REAL               RWORK( * )
  21 *       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * ),
  22 *      $                   X( LDX, * )
  23 *       ..
  24 *
  25 *
  26 *> \par Purpose:
  27 *  =============
  28 *>
  29 *> \verbatim
  30 *>
  31 *> CTRT02 computes the residual for the computed solution to a
  32 *> triangular system of linear equations  A*x = b,  A**T *x = b,
  33 *> or A**H *x = b.  Here A is a triangular matrix, A**T is the transpose
  34 *> of A, A**H is the conjugate transpose of A, and x and b are N by NRHS
  35 *> matrices.  The test ratio is the maximum over the number of right
  36 *> 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] NRHS
  76 *> \verbatim
  77 *>          NRHS is INTEGER
  78 *>          The number of right hand sides, i.e., the number of columns
  79 *>          of the matrices X and B.  NRHS >= 0.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] A
  83 *> \verbatim
  84 *>          A is COMPLEX array, dimension (LDA,N)
  85 *>          The triangular matrix A.  If UPLO = 'U', the leading n by n
  86 *>          upper triangular part of the array A contains the upper
  87 *>          triangular matrix, and the strictly lower triangular part of
  88 *>          A is not referenced.  If UPLO = 'L', the leading n by n lower
  89 *>          triangular part of the array A contains the lower triangular
  90 *>          matrix, and the strictly upper triangular part of A is not
  91 *>          referenced.  If DIAG = 'U', the diagonal elements of A are
  92 *>          also not referenced and are assumed to be 1.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] LDA
  96 *> \verbatim
  97 *>          LDA is INTEGER
  98 *>          The leading dimension of the array A.  LDA >= 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 for the system of linear
 105 *>          equations.
 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] B
 115 *> \verbatim
 116 *>          B is COMPLEX array, dimension (LDB,NRHS)
 117 *>          The right hand side vectors for the system of linear
 118 *>          equations.
 119 *> \endverbatim
 120 *>
 121 *> \param[in] LDB
 122 *> \verbatim
 123 *>          LDB is INTEGER
 124 *>          The leading dimension of the array B.  LDB >= max(1,N).
 125 *> \endverbatim
 126 *>
 127 *> \param[out] WORK
 128 *> \verbatim
 129 *>          WORK is COMPLEX array, dimension (N)
 130 *> \endverbatim
 131 *>
 132 *> \param[out] RWORK
 133 *> \verbatim
 134 *>          RWORK is REAL array, dimension (N)
 135 *> \endverbatim
 136 *>
 137 *> \param[out] RESID
 138 *> \verbatim
 139 *>          RESID is REAL
 140 *>          The maximum over the number of right hand sides of
 141 *>          norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
 142 *> \endverbatim
 143 *
 144 *  Authors:
 145 *  ========
 146 *
 147 *> \author Univ. of Tennessee
 148 *> \author Univ. of California Berkeley
 149 *> \author Univ. of Colorado Denver
 150 *> \author NAG Ltd.
 151 *
 152 *> \date November 2011
 153 *
 154 *> \ingroup complex_lin
 155 *
 156 *  =====================================================================
 157       SUBROUTINE CTRT02( UPLO, TRANS, DIAG, N, NRHS, A, LDA, X, LDX, B,
 158      $                   LDB, WORK, RWORK, RESID )
 159 *
 160 *  -- LAPACK test routine (version 3.4.0) --
 161 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 162 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 163 *     November 2011
 164 *
 165 *     .. Scalar Arguments ..
 166       CHARACTER          DIAG, TRANS, UPLO
 167       INTEGER            LDA, LDB, LDX, N, NRHS
 168       REAL               RESID
 169 *     ..
 170 *     .. Array Arguments ..
 171       REAL               RWORK( * )
 172       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * ),
 173      $                   X( LDX, * )
 174 *     ..
 175 *
 176 *  =====================================================================
 177 *
 178 *     .. Parameters ..
 179       REAL               ZERO, ONE
 180       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 181 *     ..
 182 *     .. Local Scalars ..
 183       INTEGER            J
 184       REAL               ANORM, BNORM, EPS, XNORM
 185 *     ..
 186 *     .. External Functions ..
 187       LOGICAL            LSAME
 188       REAL               CLANTR, SCASUM, SLAMCH
 189       EXTERNAL           LSAME, CLANTR, SCASUM, SLAMCH
 190 *     ..
 191 *     .. External Subroutines ..
 192       EXTERNAL           CAXPY, CCOPY, CTRMV
 193 *     ..
 194 *     .. Intrinsic Functions ..
 195       INTRINSIC          CMPLX, MAX
 196 *     ..
 197 *     .. Executable Statements ..
 198 *
 199 *     Quick exit if N = 0 or NRHS = 0
 200 *
 201       IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
 202          RESID = ZERO
 203          RETURN
 204       END IF
 205 *
 206 *     Compute the 1-norm of A or A**H.
 207 *
 208       IF( LSAME( TRANS, 'N' ) ) THEN
 209          ANORM = CLANTR( '1', UPLO, DIAG, N, N, A, LDA, RWORK )
 210       ELSE
 211          ANORM = CLANTR( 'I', UPLO, DIAG, N, N, A, LDA, RWORK )
 212       END IF
 213 *
 214 *     Exit with RESID = 1/EPS if ANORM = 0.
 215 *
 216       EPS = SLAMCH( 'Epsilon' )
 217       IF( ANORM.LE.ZERO ) THEN
 218          RESID = ONE / EPS
 219          RETURN
 220       END IF
 221 *
 222 *     Compute the maximum over the number of right hand sides of
 223 *        norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS )
 224 *
 225       RESID = ZERO
 226       DO 10 J = 1, NRHS
 227          CALL CCOPY( N, X( 1, J ), 1, WORK, 1 )
 228          CALL CTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 )
 229          CALL CAXPY( N, CMPLX( -ONE ), B( 1, J ), 1, WORK, 1 )
 230          BNORM = SCASUM( N, WORK, 1 )
 231          XNORM = SCASUM( N, X( 1, J ), 1 )
 232          IF( XNORM.LE.ZERO ) THEN
 233             RESID = ONE / EPS
 234          ELSE
 235             RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
 236          END IF
 237    10 CONTINUE
 238 *
 239       RETURN
 240 *
 241 *     End of CTRT02
 242 *
 243       END