1       SUBROUTINE DPTT02( N, NRHS, D, E, X, LDX, B, LDB, RESID )
  2 *
  3 *  -- LAPACK test routine (version 3.1) --
  4 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  5 *     November 2006
  6 *
  7 *     .. Scalar Arguments ..
  8       INTEGER            LDB, LDX, N, NRHS
  9       DOUBLE PRECISION   RESID
 10 *     ..
 11 *     .. Array Arguments ..
 12       DOUBLE PRECISION   B( LDB, * ), D( * ), E( * ), X( LDX, * )
 13 *     ..
 14 *
 15 *  Purpose
 16 *  =======
 17 *
 18 *  DPTT02 computes the residual for the solution to a symmetric
 19 *  tridiagonal system of equations:
 20 *     RESID = norm(B - A*X) / (norm(A) * norm(X) * EPS),
 21 *  where EPS is the machine epsilon.
 22 *
 23 *  Arguments
 24 *  =========
 25 *
 26 *  N       (input) INTEGTER
 27 *          The order of the matrix A.
 28 *
 29 *  NRHS    (input) INTEGER
 30 *          The number of right hand sides, i.e., the number of columns
 31 *          of the matrices B and X.  NRHS >= 0.
 32 *
 33 *  D       (input) DOUBLE PRECISION array, dimension (N)
 34 *          The n diagonal elements of the tridiagonal matrix A.
 35 *
 36 *  E       (input) DOUBLE PRECISION array, dimension (N-1)
 37 *          The (n-1) subdiagonal elements of the tridiagonal matrix A.
 38 *
 39 *  X       (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
 40 *          The n by nrhs matrix of solution vectors X.
 41 *
 42 *  LDX     (input) INTEGER
 43 *          The leading dimension of the array X.  LDX >= max(1,N).
 44 *
 45 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
 46 *          On entry, the n by nrhs matrix of right hand side vectors B.
 47 *          On exit, B is overwritten with the difference B - A*X.
 48 *
 49 *  LDB     (input) INTEGER
 50 *          The leading dimension of the array B.  LDB >= max(1,N).
 51 *
 52 *  RESID   (output) DOUBLE PRECISION
 53 *          norm(B - A*X) / (norm(A) * norm(X) * EPS)
 54 *
 55 *  =====================================================================
 56 *
 57 *     .. Parameters ..
 58       DOUBLE PRECISION   ONE, ZERO
 59       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 60 *     ..
 61 *     .. Local Scalars ..
 62       INTEGER            J
 63       DOUBLE PRECISION   ANORM, BNORM, EPS, XNORM
 64 *     ..
 65 *     .. External Functions ..
 66       DOUBLE PRECISION   DASUM, DLAMCH, DLANST
 67       EXTERNAL           DASUM, DLAMCH, DLANST
 68 *     ..
 69 *     .. Intrinsic Functions ..
 70       INTRINSIC          MAX
 71 *     ..
 72 *     .. External Subroutines ..
 73       EXTERNAL           DLAPTM
 74 *     ..
 75 *     .. Executable Statements ..
 76 *
 77 *     Quick return if possible
 78 *
 79       IF( N.LE.0 ) THEN
 80          RESID = ZERO
 81          RETURN
 82       END IF
 83 *
 84 *     Compute the 1-norm of the tridiagonal matrix A.
 85 *
 86       ANORM = DLANST( '1', N, D, E )
 87 *
 88 *     Exit with RESID = 1/EPS if ANORM = 0.
 89 *
 90       EPS = DLAMCH( 'Epsilon' )
 91       IF( ANORM.LE.ZERO ) THEN
 92          RESID = ONE / EPS
 93          RETURN
 94       END IF
 95 *
 96 *     Compute B - A*X.
 97 *
 98       CALL DLAPTM( N, NRHS, -ONE, D, E, X, LDX, ONE, B, LDB )
 99 *
100 *     Compute the maximum over the number of right hand sides of
101 *        norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
102 *
103       RESID = ZERO
104       DO 10 J = 1, NRHS
105          BNORM = DASUM( N, B( 1, J ), 1 )
106          XNORM = DASUM( N, X( 1, J ), 1 )
107          IF( XNORM.LE.ZERO ) THEN
108             RESID = ONE / EPS
109          ELSE
110             RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
111          END IF
112    10 CONTINUE
113 *
114       RETURN
115 *
116 *     End of DPTT02
117 *
118       END