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