TESTING/LIN/dtrt02.f

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