1       SUBROUTINE CSPT02( UPLO, N, NRHS, A, 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            LDB, LDX, N, NRHS
 11       REAL               RESID
 12 *     ..
 13 *     .. Array Arguments ..
 14       REAL               RWORK( * )
 15       COMPLEX            A( * ), B( LDB, * ), X( LDX, * )
 16 *     ..
 17 *
 18 *  Purpose
 19 *  =======
 20 *
 21 *  CSPT02 computes the residual in the solution of a complex symmetric
 22 *  system of linear equations  A*x = b  when packed storage is used for
 23 *  the coefficient matrix.  The ratio computed is
 24 *
 25 *     RESID = norm( B - A*X ) / ( norm(A) * norm(X) * EPS).
 26 *
 27 *  where EPS is the machine precision.
 28 *
 29 *  Arguments
 30 *  =========
 31 *
 32 *  UPLO    (input) CHARACTER*1
 33 *          Specifies whether the upper or lower triangular part of the
 34 *          complex symmetric matrix A is stored:
 35 *          = 'U':  Upper triangular
 36 *          = 'L':  Lower triangular
 37 *
 38 *  N       (input) INTEGER
 39 *          The number of rows and columns of the matrix A.  N >= 0.
 40 *
 41 *  NRHS    (input) INTEGER
 42 *          The number of columns of B, the matrix of right hand sides.
 43 *          NRHS >= 0.
 44 *
 45 *  A       (input) COMPLEX array, dimension (N*(N+1)/2)
 46 *          The original complex symmetric matrix A, stored as a packed
 47 *          triangular matrix.
 48 *
 49 *  X       (input) COMPLEX array, dimension (LDX,NRHS)
 50 *          The computed solution vectors for the system of linear
 51 *          equations.
 52 *
 53 *  LDX     (input) INTEGER
 54 *          The leading dimension of the array X.   LDX >= max(1,N).
 55 *
 56 *  B       (input/output) COMPLEX 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.  LDB >= max(1,N).
 63 *
 64 *  RWORK   (workspace) REAL array, dimension (N)
 65 *
 66 *  RESID   (output) REAL
 67 *          The maximum over the number of right hand sides of
 68 *          norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
 69 *
 70 *  =====================================================================
 71 *
 72 *     .. Parameters ..
 73       REAL               ZERO, ONE
 74       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 75       COMPLEX            CONE
 76       PARAMETER          ( CONE = ( 1.0E+00.0E+0 ) )
 77 *     ..
 78 *     .. Local Scalars ..
 79       INTEGER            J
 80       REAL               ANORM, BNORM, EPS, XNORM
 81 *     ..
 82 *     .. External Functions ..
 83       REAL               CLANSP, SCASUM, SLAMCH
 84       EXTERNAL           CLANSP, SCASUM, SLAMCH
 85 *     ..
 86 *     .. External Subroutines ..
 87       EXTERNAL           CSPMV
 88 *     ..
 89 *     .. Intrinsic Functions ..
 90       INTRINSIC          MAX
 91 *     ..
 92 *     .. Executable Statements ..
 93 *
 94 *     Quick exit if N = 0 or NRHS = 0
 95 *
 96       IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
 97          RESID = ZERO
 98          RETURN
 99       END IF
100 *
101 *     Exit with RESID = 1/EPS if ANORM = 0.
102 *
103       EPS = SLAMCH( 'Epsilon' )
104       ANORM = CLANSP( '1', UPLO, N, A, RWORK )
105       IF( ANORM.LE.ZERO ) THEN
106          RESID = ONE / EPS
107          RETURN
108       END IF
109 *
110 *     Compute  B - A*X  for the matrix of right hand sides B.
111 *
112       DO 10 J = 1, NRHS
113          CALL CSPMV( UPLO, N, -CONE, A, X( 1, J ), 1, CONE, B( 1, J ),
114      $               1 )
115    10 CONTINUE
116 *
117 *     Compute the maximum over the number of right hand sides of
118 *        norm( B - A*X ) / ( norm(A) * norm(X) * EPS ) .
119 *
120       RESID = ZERO
121       DO 20 J = 1, NRHS
122          BNORM = SCASUM( N, B( 1, J ), 1 )
123          XNORM = SCASUM( N, X( 1, J ), 1 )
124          IF( XNORM.LE.ZERO ) THEN
125             RESID = ONE / EPS
126          ELSE
127             RESID = MAX( RESID, ( ( BNORM/ANORM )/XNORM )/EPS )
128          END IF
129    20 CONTINUE
130 *
131       RETURN
132 *
133 *     End of CSPT02
134 *
135       END