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