00001 SUBROUTINE DSYT01( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC, 00002 $ RWORK, RESID ) 00003 * 00004 * -- LAPACK test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER UPLO 00010 INTEGER LDA, LDAFAC, LDC, N 00011 DOUBLE PRECISION RESID 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IPIV( * ) 00015 DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), 00016 $ RWORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * DSYT01 reconstructs a symmetric indefinite matrix A from its 00023 * block L*D*L' or U*D*U' factorization and computes the residual 00024 * norm( C - A ) / ( N * norm(A) * EPS ), 00025 * where C is the reconstructed matrix and EPS is the machine epsilon. 00026 * 00027 * Arguments 00028 * ========== 00029 * 00030 * UPLO (input) CHARACTER*1 00031 * Specifies whether the upper or lower triangular part of the 00032 * symmetric matrix A is stored: 00033 * = 'U': Upper triangular 00034 * = 'L': Lower triangular 00035 * 00036 * N (input) INTEGER 00037 * The number of rows and columns of the matrix A. N >= 0. 00038 * 00039 * A (input) DOUBLE PRECISION array, dimension (LDA,N) 00040 * The original symmetric matrix A. 00041 * 00042 * LDA (input) INTEGER 00043 * The leading dimension of the array A. LDA >= max(1,N) 00044 * 00045 * AFAC (input) DOUBLE PRECISION array, dimension (LDAFAC,N) 00046 * The factored form of the matrix A. AFAC contains the block 00047 * diagonal matrix D and the multipliers used to obtain the 00048 * factor L or U from the block L*D*L' or U*D*U' factorization 00049 * as computed by DSYTRF. 00050 * 00051 * LDAFAC (input) INTEGER 00052 * The leading dimension of the array AFAC. LDAFAC >= max(1,N). 00053 * 00054 * IPIV (input) INTEGER array, dimension (N) 00055 * The pivot indices from DSYTRF. 00056 * 00057 * C (workspace) DOUBLE PRECISION array, dimension (LDC,N) 00058 * 00059 * LDC (integer) INTEGER 00060 * The leading dimension of the array C. LDC >= max(1,N). 00061 * 00062 * RWORK (workspace) DOUBLE PRECISION array, dimension (N) 00063 * 00064 * RESID (output) DOUBLE PRECISION 00065 * If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) 00066 * If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) 00067 * 00068 * ===================================================================== 00069 * 00070 * .. Parameters .. 00071 DOUBLE PRECISION ZERO, ONE 00072 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00073 * .. 00074 * .. Local Scalars .. 00075 INTEGER I, INFO, J 00076 DOUBLE PRECISION ANORM, EPS 00077 * .. 00078 * .. External Functions .. 00079 LOGICAL LSAME 00080 DOUBLE PRECISION DLAMCH, DLANSY 00081 EXTERNAL LSAME, DLAMCH, DLANSY 00082 * .. 00083 * .. External Subroutines .. 00084 EXTERNAL DLASET, DLAVSY 00085 * .. 00086 * .. Intrinsic Functions .. 00087 INTRINSIC DBLE 00088 * .. 00089 * .. Executable Statements .. 00090 * 00091 * Quick exit if N = 0. 00092 * 00093 IF( N.LE.0 ) THEN 00094 RESID = ZERO 00095 RETURN 00096 END IF 00097 * 00098 * Determine EPS and the norm of A. 00099 * 00100 EPS = DLAMCH( 'Epsilon' ) 00101 ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK ) 00102 * 00103 * Initialize C to the identity matrix. 00104 * 00105 CALL DLASET( 'Full', N, N, ZERO, ONE, C, LDC ) 00106 * 00107 * Call DLAVSY to form the product D * U' (or D * L' ). 00108 * 00109 CALL DLAVSY( UPLO, 'Transpose', 'Non-unit', N, N, AFAC, LDAFAC, 00110 $ IPIV, C, LDC, INFO ) 00111 * 00112 * Call DLAVSY again to multiply by U (or L ). 00113 * 00114 CALL DLAVSY( UPLO, 'No transpose', 'Unit', N, N, AFAC, LDAFAC, 00115 $ IPIV, C, LDC, INFO ) 00116 * 00117 * Compute the difference C - A . 00118 * 00119 IF( LSAME( UPLO, 'U' ) ) THEN 00120 DO 20 J = 1, N 00121 DO 10 I = 1, J 00122 C( I, J ) = C( I, J ) - A( I, J ) 00123 10 CONTINUE 00124 20 CONTINUE 00125 ELSE 00126 DO 40 J = 1, N 00127 DO 30 I = J, N 00128 C( I, J ) = C( I, J ) - A( I, J ) 00129 30 CONTINUE 00130 40 CONTINUE 00131 END IF 00132 * 00133 * Compute norm( C - A ) / ( N * norm(A) * EPS ) 00134 * 00135 RESID = DLANSY( '1', UPLO, N, C, LDC, RWORK ) 00136 * 00137 IF( ANORM.LE.ZERO ) THEN 00138 IF( RESID.NE.ZERO ) 00139 $ RESID = ONE / EPS 00140 ELSE 00141 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS 00142 END IF 00143 * 00144 RETURN 00145 * 00146 * End of DSYT01 00147 * 00148 END
1.7.2