![]() |
LAPACK 3.3.1
Linear Algebra PACKage
|
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