1*> \brief \b ZSPT01 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 ZSPT01( UPLO, N, A, AFAC, IPIV, C, LDC, RWORK, RESID ) 12* 13* .. Scalar Arguments .. 14* CHARACTER UPLO 15* INTEGER LDC, N 16* DOUBLE PRECISION RESID 17* .. 18* .. Array Arguments .. 19* INTEGER IPIV( * ) 20* DOUBLE PRECISION RWORK( * ) 21* COMPLEX*16 A( * ), AFAC( * ), C( LDC, * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> ZSPT01 reconstructs a symmetric indefinite packed matrix A from its 31*> diagonal pivoting factorization A = U*D*U' or A = L*D*L' and computes 32*> the residual 33*> norm( C - A ) / ( N * norm(A) * EPS ), 34*> where C is the reconstructed matrix and EPS is the machine epsilon. 35*> \endverbatim 36* 37* Arguments: 38* ========== 39* 40*> \param[in] UPLO 41*> \verbatim 42*> UPLO is CHARACTER*1 43*> Specifies whether the upper or lower triangular part of the 44*> Hermitian matrix A is stored: 45*> = 'U': Upper triangular 46*> = 'L': Lower triangular 47*> \endverbatim 48*> 49*> \param[in] N 50*> \verbatim 51*> N is INTEGER 52*> The order of the matrix A. N >= 0. 53*> \endverbatim 54*> 55*> \param[in] A 56*> \verbatim 57*> A is COMPLEX*16 array, dimension (N*(N+1)/2) 58*> The original symmetric matrix A, stored as a packed 59*> triangular matrix. 60*> \endverbatim 61*> 62*> \param[in] AFAC 63*> \verbatim 64*> AFAC is COMPLEX*16 array, dimension (N*(N+1)/2) 65*> The factored form of the matrix A, stored as a packed 66*> triangular matrix. AFAC contains the block diagonal matrix D 67*> and the multipliers used to obtain the factor L or U from the 68*> L*D*L' or U*D*U' factorization as computed by ZSPTRF. 69*> \endverbatim 70*> 71*> \param[in] IPIV 72*> \verbatim 73*> IPIV is INTEGER array, dimension (N) 74*> The pivot indices from ZSPTRF. 75*> \endverbatim 76*> 77*> \param[out] C 78*> \verbatim 79*> C is COMPLEX*16 array, dimension (LDC,N) 80*> \endverbatim 81*> 82*> \param[in] LDC 83*> \verbatim 84*> LDC is INTEGER 85*> The leading dimension of the array C. LDC >= max(1,N). 86*> \endverbatim 87*> 88*> \param[out] RWORK 89*> \verbatim 90*> RWORK is DOUBLE PRECISION array, dimension (N) 91*> \endverbatim 92*> 93*> \param[out] RESID 94*> \verbatim 95*> RESID is DOUBLE PRECISION 96*> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) 97*> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) 98*> \endverbatim 99* 100* Authors: 101* ======== 102* 103*> \author Univ. of Tennessee 104*> \author Univ. of California Berkeley 105*> \author Univ. of Colorado Denver 106*> \author NAG Ltd. 107* 108*> \date December 2016 109* 110*> \ingroup complex16_lin 111* 112* ===================================================================== 113 SUBROUTINE ZSPT01( UPLO, N, A, AFAC, IPIV, C, LDC, RWORK, RESID ) 114* 115* -- LAPACK test routine (version 3.7.0) -- 116* -- LAPACK is a software package provided by Univ. of Tennessee, -- 117* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 118* December 2016 119* 120* .. Scalar Arguments .. 121 CHARACTER UPLO 122 INTEGER LDC, N 123 DOUBLE PRECISION RESID 124* .. 125* .. Array Arguments .. 126 INTEGER IPIV( * ) 127 DOUBLE PRECISION RWORK( * ) 128 COMPLEX*16 A( * ), AFAC( * ), C( LDC, * ) 129* .. 130* 131* ===================================================================== 132* 133* .. Parameters .. 134 DOUBLE PRECISION ZERO, ONE 135 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 136 COMPLEX*16 CZERO, CONE 137 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 138 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 139* .. 140* .. Local Scalars .. 141 INTEGER I, INFO, J, JC 142 DOUBLE PRECISION ANORM, EPS 143* .. 144* .. External Functions .. 145 LOGICAL LSAME 146 DOUBLE PRECISION DLAMCH, ZLANSP, ZLANSY 147 EXTERNAL LSAME, DLAMCH, ZLANSP, ZLANSY 148* .. 149* .. External Subroutines .. 150 EXTERNAL ZLASET, ZLAVSP 151* .. 152* .. Intrinsic Functions .. 153 INTRINSIC DBLE 154* .. 155* .. Executable Statements .. 156* 157* Quick exit if N = 0. 158* 159 IF( N.LE.0 ) THEN 160 RESID = ZERO 161 RETURN 162 END IF 163* 164* Determine EPS and the norm of A. 165* 166 EPS = DLAMCH( 'Epsilon' ) 167 ANORM = ZLANSP( '1', UPLO, N, A, RWORK ) 168* 169* Initialize C to the identity matrix. 170* 171 CALL ZLASET( 'Full', N, N, CZERO, CONE, C, LDC ) 172* 173* Call ZLAVSP to form the product D * U' (or D * L' ). 174* 175 CALL ZLAVSP( UPLO, 'Transpose', 'Non-unit', N, N, AFAC, IPIV, C, 176 $ LDC, INFO ) 177* 178* Call ZLAVSP again to multiply by U ( or L ). 179* 180 CALL ZLAVSP( UPLO, 'No transpose', 'Unit', N, N, AFAC, IPIV, C, 181 $ LDC, INFO ) 182* 183* Compute the difference C - A . 184* 185 IF( LSAME( UPLO, 'U' ) ) THEN 186 JC = 0 187 DO 20 J = 1, N 188 DO 10 I = 1, J 189 C( I, J ) = C( I, J ) - A( JC+I ) 190 10 CONTINUE 191 JC = JC + J 192 20 CONTINUE 193 ELSE 194 JC = 1 195 DO 40 J = 1, N 196 DO 30 I = J, N 197 C( I, J ) = C( I, J ) - A( JC+I-J ) 198 30 CONTINUE 199 JC = JC + N - J + 1 200 40 CONTINUE 201 END IF 202* 203* Compute norm( C - A ) / ( N * norm(A) * EPS ) 204* 205 RESID = ZLANSY( '1', UPLO, N, C, LDC, RWORK ) 206* 207 IF( ANORM.LE.ZERO ) THEN 208 IF( RESID.NE.ZERO ) 209 $ RESID = ONE / EPS 210 ELSE 211 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS 212 END IF 213* 214 RETURN 215* 216* End of ZSPT01 217* 218 END 219