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