d8/d35/cptt01_8f_source.html

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