TESTING/LIN/dpot02.f

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