da/d7e/dptt01_8f_source.html

00001       SUBROUTINE DPTT01( N, D, E, DF, EF, WORK, 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            N
00009       DOUBLE PRECISION   RESID
00010 *     ..
00011 *     .. Array Arguments ..
00012       DOUBLE PRECISION   D( * ), DF( * ), E( * ), EF( * ), WORK( * )
00013 *     ..
00014 *
00015 *  Purpose
00016 *  =======
00017 *
00018 *  DPTT01 reconstructs a tridiagonal matrix A from its L*D*L'
00019 *  factorization and computes the residual
00020 *     norm(L*D*L' - A) / ( n * norm(A) * 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 *  D       (input) DOUBLE PRECISION array, dimension (N)
00030 *          The n diagonal elements of the tridiagonal matrix A.
00031 *
00032 *  E       (input) DOUBLE PRECISION array, dimension (N-1)
00033 *          The (n-1) subdiagonal elements of the tridiagonal matrix A.
00034 *
00035 *  DF      (input) DOUBLE PRECISION array, dimension (N)
00036 *          The n diagonal elements of the factor L from the L*D*L'
00037 *          factorization of A.
00038 *
00039 *  EF      (input) DOUBLE PRECISION array, dimension (N-1)
00040 *          The (n-1) subdiagonal elements of the factor L from the
00041 *          L*D*L' factorization of A.
00042 *
00043 *  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)
00044 *
00045 *  RESID   (output) DOUBLE PRECISION
00046 *          norm(L*D*L' - A) / (n * norm(A) * EPS)
00047 *
00048 *  =====================================================================
00049 *
00050 *     .. Parameters ..
00051       DOUBLE PRECISION   ONE, ZERO
00052       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00053 *     ..
00054 *     .. Local Scalars ..
00055       INTEGER            I
00056       DOUBLE PRECISION   ANORM, DE, EPS
00057 *     ..
00058 *     .. External Functions ..
00059       DOUBLE PRECISION   DLAMCH
00060       EXTERNAL           DLAMCH
00061 *     ..
00062 *     .. Intrinsic Functions ..
00063       INTRINSIC          ABS, DBLE, MAX
00064 *     ..
00065 *     .. Executable Statements ..
00066 *
00067 *     Quick return if possible
00068 *
00069       IF( N.LE.0 ) THEN
00070          RESID = ZERO
00071          RETURN
00072       END IF
00073 *
00074       EPS = DLAMCH( 'Epsilon' )
00075 *
00076 *     Construct the difference L*D*L' - A.
00077 *
00078       WORK( 1 ) = DF( 1 ) - D( 1 )
00079       DO 10 I = 1, N - 1
00080          DE = DF( I )*EF( I )
00081          WORK( N+I ) = DE - E( I )
00082          WORK( 1+I ) = DE*EF( I ) + DF( I+1 ) - D( I+1 )
00083    10 CONTINUE
00084 *
00085 *     Compute the 1-norms of the tridiagonal matrices A and WORK.
00086 *
00087       IF( N.EQ.1 ) THEN
00088          ANORM = D( 1 )
00089          RESID = ABS( WORK( 1 ) )
00090       ELSE
00091          ANORM = MAX( D( 1 )+ABS( E( 1 ) ), D( N )+ABS( E( N-1 ) ) )
00092          RESID = MAX( ABS( WORK( 1 ) )+ABS( WORK( N+1 ) ),
00093      $           ABS( WORK( N ) )+ABS( WORK( 2*N-1 ) ) )
00094          DO 20 I = 2, N - 1
00095             ANORM = MAX( ANORM, D( I )+ABS( E( I ) )+ABS( E( I-1 ) ) )
00096             RESID = MAX( RESID, ABS( WORK( I ) )+ABS( WORK( N+I-1 ) )+
00097      $              ABS( WORK( N+I ) ) )
00098    20    CONTINUE
00099       END IF
00100 *
00101 *     Compute norm(L*D*L' - A) / (n * norm(A) * EPS)
00102 *
00103       IF( ANORM.LE.ZERO ) THEN
00104          IF( RESID.NE.ZERO )
00105      $      RESID = ONE / EPS
00106       ELSE
00107          RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
00108       END IF
00109 *
00110       RETURN
00111 *
00112 *     End of DPTT01
00113 *
00114       END