1*> \brief \b DSYT01_3 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE DSYT01_3( UPLO, N, A, LDA, AFAC, LDAFAC, E, IPIV, C, 12* LDC, RWORK, RESID ) 13* 14* .. Scalar Arguments .. 15* CHARACTER UPLO 16* INTEGER LDA, LDAFAC, LDC, N 17* DOUBLE PRECISION RESID 18* .. 19* .. Array Arguments .. 20* INTEGER IPIV( * ) 21* DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), 22* $ E( * ), RWORK( * ) 23* .. 24* 25* 26*> \par Purpose: 27* ============= 28*> 29*> \verbatim 30*> 31*> DSYT01_3 reconstructs a symmetric indefinite matrix A from its 32*> block L*D*L' or U*D*U' factorization computed by DSYTRF_RK 33*> (or DSYTRF_BK) and computes the residual 34*> norm( C - A ) / ( N * norm(A) * EPS ), 35*> where C is the reconstructed matrix and EPS is the machine epsilon. 36*> \endverbatim 37* 38* Arguments: 39* ========== 40* 41*> \param[in] UPLO 42*> \verbatim 43*> UPLO is CHARACTER*1 44*> Specifies whether the upper or lower triangular part of the 45*> symmetric matrix A is stored: 46*> = 'U': Upper triangular 47*> = 'L': Lower triangular 48*> \endverbatim 49*> 50*> \param[in] N 51*> \verbatim 52*> N is INTEGER 53*> The number of rows and columns of the matrix A. N >= 0. 54*> \endverbatim 55*> 56*> \param[in] A 57*> \verbatim 58*> A is DOUBLE PRECISION array, dimension (LDA,N) 59*> The original symmetric matrix A. 60*> \endverbatim 61*> 62*> \param[in] LDA 63*> \verbatim 64*> LDA is INTEGER 65*> The leading dimension of the array A. LDA >= max(1,N) 66*> \endverbatim 67*> 68*> \param[in] AFAC 69*> \verbatim 70*> AFAC is DOUBLE PRECISION array, dimension (LDAFAC,N) 71*> Diagonal of the block diagonal matrix D and factors U or L 72*> as computed by DSYTRF_RK and DSYTRF_BK: 73*> a) ONLY diagonal elements of the symmetric block diagonal 74*> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k); 75*> (superdiagonal (or subdiagonal) elements of D 76*> should be provided on entry in array E), and 77*> b) If UPLO = 'U': factor U in the superdiagonal part of A. 78*> If UPLO = 'L': factor L in the subdiagonal part of A. 79*> \endverbatim 80*> 81*> \param[in] LDAFAC 82*> \verbatim 83*> LDAFAC is INTEGER 84*> The leading dimension of the array AFAC. 85*> LDAFAC >= max(1,N). 86*> \endverbatim 87*> 88*> \param[in] E 89*> \verbatim 90*> E is DOUBLE PRECISION array, dimension (N) 91*> On entry, contains the superdiagonal (or subdiagonal) 92*> elements of the symmetric block diagonal matrix D 93*> with 1-by-1 or 2-by-2 diagonal blocks, where 94*> If UPLO = 'U': E(i) = D(i-1,i),i=2:N, E(1) not referenced; 95*> If UPLO = 'L': E(i) = D(i+1,i),i=1:N-1, E(N) not referenced. 96*> \endverbatim 97*> 98*> \param[in] IPIV 99*> \verbatim 100*> IPIV is INTEGER array, dimension (N) 101*> The pivot indices from DSYTRF_RK (or DSYTRF_BK). 102*> \endverbatim 103*> 104*> \param[out] C 105*> \verbatim 106*> C is DOUBLE PRECISION array, dimension (LDC,N) 107*> \endverbatim 108*> 109*> \param[in] LDC 110*> \verbatim 111*> LDC is INTEGER 112*> The leading dimension of the array C. LDC >= max(1,N). 113*> \endverbatim 114*> 115*> \param[out] RWORK 116*> \verbatim 117*> RWORK is DOUBLE PRECISION array, dimension (N) 118*> \endverbatim 119*> 120*> \param[out] RESID 121*> \verbatim 122*> RESID is DOUBLE PRECISION 123*> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) 124*> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) 125*> \endverbatim 126* 127* Authors: 128* ======== 129* 130*> \author Univ. of Tennessee 131*> \author Univ. of California Berkeley 132*> \author Univ. of Colorado Denver 133*> \author NAG Ltd. 134* 135*> \date June 2017 136* 137*> \ingroup double_lin 138* 139* ===================================================================== 140 SUBROUTINE DSYT01_3( UPLO, N, A, LDA, AFAC, LDAFAC, E, IPIV, C, 141 $ LDC, RWORK, RESID ) 142* 143* -- LAPACK test routine (version 3.7.1) -- 144* -- LAPACK is a software package provided by Univ. of Tennessee, -- 145* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 146* June 2017 147* 148* .. Scalar Arguments .. 149 CHARACTER UPLO 150 INTEGER LDA, LDAFAC, LDC, N 151 DOUBLE PRECISION RESID 152* .. 153* .. Array Arguments .. 154 INTEGER IPIV( * ) 155 DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), 156 $ E( * ), RWORK( * ) 157* .. 158* 159* ===================================================================== 160* 161* .. Parameters .. 162 DOUBLE PRECISION ZERO, ONE 163 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 164* .. 165* .. Local Scalars .. 166 INTEGER I, INFO, J 167 DOUBLE PRECISION ANORM, EPS 168* .. 169* .. External Functions .. 170 LOGICAL LSAME 171 DOUBLE PRECISION DLAMCH, DLANSY 172 EXTERNAL LSAME, DLAMCH, DLANSY 173* .. 174* .. External Subroutines .. 175 EXTERNAL DLASET, DLAVSY_ROOK, DSYCONVF_ROOK 176* .. 177* .. Intrinsic Functions .. 178 INTRINSIC DBLE 179* .. 180* .. Executable Statements .. 181* 182* Quick exit if N = 0. 183* 184 IF( N.LE.0 ) THEN 185 RESID = ZERO 186 RETURN 187 END IF 188* 189* a) Revert to multiplyers of L 190* 191 CALL DSYCONVF_ROOK( UPLO, 'R', N, AFAC, LDAFAC, E, IPIV, INFO ) 192* 193* 1) Determine EPS and the norm of A. 194* 195 EPS = DLAMCH( 'Epsilon' ) 196 ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK ) 197* 198* 2) Initialize C to the identity matrix. 199* 200 CALL DLASET( 'Full', N, N, ZERO, ONE, C, LDC ) 201* 202* 3) Call DLAVSY_ROOK to form the product D * U' (or D * L' ). 203* 204 CALL DLAVSY_ROOK( UPLO, 'Transpose', 'Non-unit', N, N, AFAC, 205 $ LDAFAC, IPIV, C, LDC, INFO ) 206* 207* 4) Call DLAVSY_ROOK again to multiply by U (or L ). 208* 209 CALL DLAVSY_ROOK( UPLO, 'No transpose', 'Unit', N, N, AFAC, 210 $ LDAFAC, IPIV, C, LDC, INFO ) 211* 212* 5) Compute the difference C - A. 213* 214 IF( LSAME( UPLO, 'U' ) ) THEN 215 DO J = 1, N 216 DO I = 1, J 217 C( I, J ) = C( I, J ) - A( I, J ) 218 END DO 219 END DO 220 ELSE 221 DO J = 1, N 222 DO I = J, N 223 C( I, J ) = C( I, J ) - A( I, J ) 224 END DO 225 END DO 226 END IF 227* 228* 6) Compute norm( C - A ) / ( N * norm(A) * EPS ) 229* 230 RESID = DLANSY( '1', UPLO, N, C, LDC, RWORK ) 231* 232 IF( ANORM.LE.ZERO ) THEN 233 IF( RESID.NE.ZERO ) 234 $ RESID = ONE / EPS 235 ELSE 236 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS 237 END IF 238 239* 240* b) Convert to factor of L (or U) 241* 242 CALL DSYCONVF_ROOK( UPLO, 'C', N, AFAC, LDAFAC, E, IPIV, INFO ) 243* 244 RETURN 245* 246* End of DSYT01_3 247* 248 END 249