d7/dbd/dptt02_8f_source.html

00001       SUBROUTINE DPTT02( N, NRHS, D, E, X, LDX, B, LDB, RESID )
00002 *
00003 *  -- LAPACK test routine (version 3.1) --
00004 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00005 *     November 2006
00006 *
00007 *     .. Scalar Arguments ..
00008       INTEGER            LDB, LDX, N, NRHS
00009       DOUBLE PRECISION   RESID
00010 *     ..
00011 *     .. Array Arguments ..
00012       DOUBLE PRECISION   B( LDB, * ), D( * ), E( * ), X( LDX, * )
00013 *     ..
00014 *
00015 *  Purpose
00016 *  =======
00017 *
00018 *  DPTT02 computes the residual for the solution to a symmetric
00019 *  tridiagonal system of equations:
00020 *     RESID = norm(B - A*X) / (norm(A) * norm(X) * EPS),
00021 *  where EPS is the machine epsilon.
00022 *
00023 *  Arguments
00024 *  =========
00025 *
00026 *  N       (input) INTEGTER
00027 *          The order of the matrix A.
00028 *
00029 *  NRHS    (input) INTEGER
00030 *          The number of right hand sides, i.e., the number of columns
00031 *          of the matrices B and X.  NRHS >= 0.
00032 *
00033 *  D       (input) DOUBLE PRECISION array, dimension (N)
00034 *          The n diagonal elements of the tridiagonal matrix A.
00035 *
00036 *  E       (input) DOUBLE PRECISION array, dimension (N-1)
00037 *          The (n-1) subdiagonal elements of the tridiagonal matrix A.
00038 *
00039 *  X       (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
00040 *          The n by nrhs matrix of solution vectors X.
00041 *
00042 *  LDX     (input) INTEGER
00043 *          The leading dimension of the array X.  LDX >= max(1,N).
00044 *
00045 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
00046 *          On entry, the n by nrhs matrix of right hand side vectors B.
00047 *          On exit, B is overwritten with the difference B - A*X.
00048 *
00049 *  LDB     (input) INTEGER
00050 *          The leading dimension of the array B.  LDB >= max(1,N).
00051 *
00052 *  RESID   (output) DOUBLE PRECISION
00053 *          norm(B - A*X) / (norm(A) * norm(X) * EPS)
00054 *
00055 *  =====================================================================
00056 *
00057 *     .. Parameters ..
00058       DOUBLE PRECISION   ONE, ZERO
00059       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00060 *     ..
00061 *     .. Local Scalars ..
00062       INTEGER            J
00063       DOUBLE PRECISION   ANORM, BNORM, EPS, XNORM
00064 *     ..
00065 *     .. External Functions ..
00066       DOUBLE PRECISION   DASUM, DLAMCH, DLANST
00067       EXTERNAL           DASUM, DLAMCH, DLANST
00068 *     ..
00069 *     .. Intrinsic Functions ..
00070       INTRINSIC          MAX
00071 *     ..
00072 *     .. External Subroutines ..
00073       EXTERNAL           DLAPTM
00074 *     ..
00075 *     .. Executable Statements ..
00076 *
00077 *     Quick return if possible
00078 *
00079       IF( N.LE.0 ) THEN
00080          RESID = ZERO
00081          RETURN
00082       END IF
00083 *
00084 *     Compute the 1-norm of the tridiagonal matrix A.
00085 *
00086       ANORM = DLANST( '1', N, D, E )
00087 *
00088 *     Exit with RESID = 1/EPS if ANORM = 0.
00089 *
00090       EPS = DLAMCH( 'Epsilon' )
00091       IF( ANORM.LE.ZERO ) THEN
00092          RESID = ONE / EPS
00093          RETURN
00094       END IF
00095 *
00096 *     Compute B - A*X.
00097 *
00098       CALL DLAPTM( N, NRHS, -ONE, D, E, X, LDX, ONE, B, LDB )
00099 *
00100 *     Compute the maximum over the number of right hand sides of
00101 *        norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
00102 *
00103       RESID = ZERO
00104       DO 10 J = 1, NRHS
00105          BNORM = DASUM( N, B( 1, J ), 1 )
00106          XNORM = DASUM( N, X( 1, J ), 1 )
00107          IF( XNORM.LE.ZERO ) THEN
00108             RESID = ONE / EPS
00109          ELSE
00110             RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
00111          END IF
00112    10 CONTINUE
00113 *
00114       RETURN
00115 *
00116 *     End of DPTT02
00117 *
00118       END