1       SUBROUTINE DPOT06( UPLO, N, NRHS, A, LDA, X, LDX, B, LDB,
  2      $                   RWORK, RESID )
  3 *
  4 *  -- LAPACK test routine (version 3.1.2) --
  5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  6 *     April 2007
  7 *
  8 *     .. Scalar Arguments ..
  9       CHARACTER          UPLO
 10       INTEGER            LDA, LDB, LDX, N, NRHS
 11       DOUBLE PRECISION   RESID
 12 *     ..
 13 *     .. Array Arguments ..
 14       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), RWORK( * ),
 15      $                   X( LDX, * )
 16 *     ..
 17 *
 18 *  Purpose
 19 *  =======
 20 *
 21 *  DPOT06 computes the residual for a solution of a system of linear
 22 *  equations  A*x = b :
 23 *     RESID = norm(B - A*X,inf) / ( norm(A,inf) * norm(X,inf) * EPS ),
 24 *  where EPS is the machine epsilon.
 25 *
 26 *  Arguments
 27 *  =========
 28 *
 29 *  UPLO    (input) CHARACTER*1
 30 *          Specifies whether the upper or lower triangular part of the
 31 *          symmetric matrix A is stored:
 32 *          = 'U':  Upper triangular
 33 *          = 'L':  Lower triangular
 34 *
 35 *  N       (input) INTEGER
 36 *          The number of rows and columns of the matrix A.  N >= 0.
 37 *
 38 *  NRHS    (input) INTEGER
 39 *          The number of columns of B, the matrix of right hand sides.
 40 *          NRHS >= 0.
 41 *
 42 *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
 43 *          The original M x N matrix A.
 44 *
 45 *  LDA     (input) INTEGER
 46 *          The leading dimension of the array A.  LDA >= max(1,N).
 47 *
 48 *  X       (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
 49 *          The computed solution vectors for the system of linear
 50 *          equations.
 51 *
 52 *  LDX     (input) INTEGER
 53 *          The leading dimension of the array X.  If TRANS = 'N',
 54 *          LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,N).
 55 *
 56 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
 57 *          On entry, the right hand side vectors for the system of
 58 *          linear equations.
 59 *          On exit, B is overwritten with the difference B - A*X.
 60 *
 61 *  LDB     (input) INTEGER
 62 *          The leading dimension of the array B.  IF TRANS = 'N',
 63 *          LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N).
 64 *
 65 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
 66 *
 67 *  RESID   (output) DOUBLE PRECISION
 68 *          The maximum over the number of right hand sides of
 69 *          norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
 70 *
 71 *  =====================================================================
 72 *
 73 *     .. Parameters ..
 74       DOUBLE PRECISION   ZERO, ONE, NEGONE
 75       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 76       PARAMETER          ( NEGONE = -1.0D+0 )
 77 *     ..
 78 *     .. Local Scalars ..
 79       INTEGER            IFAIL, J
 80       DOUBLE PRECISION   ANORM, BNORM, EPS, XNORM
 81 *     ..
 82 *     .. External Functions ..
 83       LOGICAL            LSAME
 84       INTEGER            IDAMAX
 85       DOUBLE PRECISION   DLAMCH, DLANSY
 86       EXTERNAL           LSAME, IDAMAX, DLAMCH, DLANSY
 87 *     ..
 88 *     .. External Subroutines ..
 89       EXTERNAL           DSYMM
 90 *     ..
 91 *     .. Intrinsic Functions ..
 92       INTRINSIC          MAXABS
 93 *     ..
 94 *     .. Executable Statements ..
 95 *
 96 *     Quick exit if N = 0 or NRHS = 0
 97 *
 98       IF( N.LE.0 .OR. NRHS.EQ.0 ) THEN
 99          RESID = ZERO
100          RETURN
101       END IF
102 *
103 *     Exit with RESID = 1/EPS if ANORM = 0.
104 *
105       EPS = DLAMCH( 'Epsilon' )
106       ANORM = DLANSY( 'I', UPLO, N, A, LDA, RWORK )
107       IF( ANORM.LE.ZERO ) THEN
108          RESID = ONE / EPS
109          RETURN
110       END IF
111 *
112 *     Compute  B - A*X  and store in B.
113       IFAIL=0
114 *
115       CALL DSYMM( 'Left', UPLO, N, NRHS, NEGONE, A, LDA, X,
116      $            LDX, ONE, B, LDB )
117 *
118 *     Compute the maximum over the number of right hand sides of
119 *        norm(B - A*X) / ( norm(A) * norm(X) * EPS ) .
120 *
121       RESID = ZERO
122       DO 10 J = 1, NRHS
123          BNORM = ABS(B(IDAMAX( N, B( 1, J ), 1 ),J))
124          XNORM = ABS(X(IDAMAX( N, X( 1, J ), 1 ),J))
125          IF( XNORM.LE.ZERO ) THEN
126             RESID = ONE / EPS
127          ELSE
128             RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
129          END IF
130    10 CONTINUE
131 *
132       RETURN
133 *
134 *     End of DPOT06
135 *
136       END